Supplemental materials for: Rolek, B.W., McClure, CJW, Dunn, L., Curti, M., … Ridgway’s Hawk IPM and PVA

Contact information:

Metadata, data, and scripts used in analyses can be found at https://github.com/The-Peregrine-Fund/XXXXX.

The full workflow below is visible as a html website at: https://the-peregrine-fund.github.io/XXXXX/.

A permanent archive and DOI is available at: https://zenodo.org/doi/XXXXX


1. Code for IPM and PVA

1.1 Submodels without Integration

library('nimble')
library('parallel')
library ('coda')
load("/bsuscratch/brianrolek/riha_ipm/data.rdata")

#**********************
#* Parameter descriptions
#**********************
# Mark-resight-recovery data
#   Observations (po) = y  
#     1 seen first-year (age=0, just before 1st b-day)
#     2 seen nonbreeder
#     3 seen breeder
#     4 not seen
#   States (ps)
#     1 alive first-year
#     2 alive nonbreeder
#     3 alive breeder
#     4 dead
#   Groups
#     1 wild-born
#     2 translocated and hacked
###################################################
# PARAMETERS
#   phiFY: survival probability first year 
#   phiA: survival probability nonbreeders
#   phiB: survival probability breeders
#   psiFYB: recruitment probability from first-year to breeder
#   psiAB: recruitment probability from nonbreeders to breeder
#   psiBA: recruitment probability from breeder to nonbreeders 
#   pFY: resight probability first-years
#   pA: resight probability nonbreeders
#   pB: resight probability breeders
#   lmu.prod: mean productivity (males and females) per territory on the log scale
#   sds: standard deviations for multivariate normal random effects for time and site
#   R: correlation coefficients for multivariate normal random effects for time and site
#   lambda: population growth rate (derived)
#   extinct: binary indicator of extirpation at a site
#   gamma: coefficient of effect from nest treatments
#   betas: coefficient of effect from translocations
#   deltas: coefficient of effect from survey effort
#   mus: overall means for survival, recruitment, and detections
#   r and rr: "r" parameter for negative binomial distribution
#             also described as omega in manuscript

#**********************
#* Model code
#**********************
mycode <- nimbleCode(
  {
    ###################################################
    # Priors and constraints
    ###################################################
    # Survival, recruitment, and detection can be correlated
    for (k in 1:8){ 
      betas[k] ~ dunif(-20, 20)  # prior for coefficients
    } # k
    for (j in 1:8){ 
      for (s in 1:nsite){
        lmus[j,s] <- logit(mus[j,s])
        mus[j,s] ~ dbeta(1,1) # prior for means
      }}     # m population #s sex #h hacked 
    
    # Temporal random effects and correlations among all sites, synchrony
    for (j in 1:p){ sds[j] ~ dexp(1) }# prior for temporal variation estimated using the multivariate normal distribution
    R[1:p,1:p] <- t(Ustar[1:p,1:p]) %*% Ustar[1:p,1:p] # calculate rhos, correlation coefficients
    Ustar[1:p,1:p] ~ dlkj_corr_cholesky(eta=1.1, p=p) # Ustar is the Cholesky decomposition of the correlation matrix
    U[1:p,1:p] <- uppertri_mult_diag(Ustar[1:p, 1:p], sds[1:p])
    # Multivariate normal for temporal variance
    for (t in 1:nyr){ 
      eps[1:p,t] ~ dmnorm(mu.zeroes[1:p],
                          cholesky = U[1:p, 1:p], prec_param = 0)
    }
    
    # Temporal random effects and correlations between sites
    for (jj in 1:p2){ sds2[jj] ~ dexp(1) }# prior for temporal variation estimated using the multivariate normal distribution
    R2[1:p2,1:p2] <- t(Ustar2[1:p2,1:p2]) %*% Ustar2[1:p2,1:p2] # calculate rhos, correlation coefficients
    Ustar2[1:p2,1:p2] ~ dlkj_corr_cholesky(eta=1.1, p=p2) # Ustar is the Cholesky decomposition of the correlation matrix
    U2[1:p2,1:p2] <- uppertri_mult_diag(Ustar2[1:p2, 1:p2], sds2[1:p2])
    # multivariate normal for temporal and site variance
    for (t in 1:nyr){ 
      for (s in 1:nsite){
        eta[1:p2,s,t] ~ dmnorm(mu.zeroes2[1:p2],
                               cholesky = U2[1:p2, 1:p2], prec_param = 0)
      } } # s t 

    ###############################
    # Likelihood for productivity
    ###############################
    # Priors for productivity
    for (s in 1:nsite){
      lmu.f[s] ~ dnorm(0, sd=5)
    } # s
    gamma ~ dunif(-20, 20)
    rr ~ dexp(0.05)
    
    # Productivity likelihood      
    for (k in 1:nnest){
      f[k] ~ dnegbin(ppp[k], rr)
      ppp[k] <- rr/(rr+mu.f[k])
      log(mu.f[k]) <- lmu.f[site.nest[k]] +  
        gamma*treat.nest[k] + 
        eps[9, year.nest[k] ] + 
        eta[9, site.nest[k], year.nest[k] ] 
    } # k
    # Derive yearly brood size for population model
    # Need to reorder because nimble doesn't 
    # handle nonconsecutive indices
    # yrind.nest is a matrix of indices for each site
    for (t in 1:nyr){
      for (s in 1:nsite){
        for (xxx in 1:nest.end[t,s]){
          fecmat[t,s,xxx] <- mu.f[ yrind.nest[xxx,t,s] ]
        } # xxx
        mn.f[t,s] <- mean( fecmat[t,s,1:nest.end[t,s]] )
      }} # s t
    
    # GOF for number of fledglings
    for (k in 1:nnest){
      f.obs[k] <- f[k] # observed counts
      f.exp[k] <- mu.f[k] # expected counts adult breeder
      f.rep[k] ~ dnegbin(ppp[k], rr) # expected counts
      f.dssm.obs[k] <- abs( ( f.obs[k] - f.exp[k] ) / (f.obs[k]+0.001) )
      f.dssm.rep[k] <- abs( ( f.rep[k] - f.exp[k] ) / (f.rep[k]+0.001) )
    } # k
    f.dmape.obs <- sum(f.dssm.obs[1:nnest])
    f.dmape.rep <- sum(f.dssm.rep[1:nnest])
    f.tvm.obs <- sd(brood[1:nnest])^2/mean(brood[1:nnest])
    f.tvm.rep <- sd(f.rep[1:nnest])^2/mean(f.rep[1:nnest])
    
    ################################
    # Likelihood for counts
    ################################
# Omitted. See IPM or PVA.
    
    ################################
    # Likelihood for survival
    ################################ 
    # Calculate yearly averages for sites for integration
    for (t in 1:nyr){
      for (s in 1:nsite){
        for (xxxx in 1:surv.end[t,s]){
          # Need to reorder because nimble doesn't 
          # handle nonconsecutive indices
          # yrind.surv is a matrix of indices for each site
          phiFY2[ t, s, xxxx] <- phiFY[ yrind.surv[xxxx,t,s], t]
          phiA2[ t, s, xxxx] <- phiA[ yrind.surv[xxxx,t,s], t]
          phiB2[ t, s, xxxx] <- phiB[ yrind.surv[xxxx,t,s], t]
          psiFYB2[ t, s, xxxx] <- psiFYB[ yrind.surv[xxxx,t,s], t]
          psiAB2[ t, s, xxxx] <- psiAB[ yrind.surv[xxxx,t,s], t]
          psiBA2[ t, s, xxxx] <- psiBA[ yrind.surv[xxxx,t,s], t]
          pA2[ t, s, xxxx] <- pA[ yrind.surv[xxxx,t,s], t]
          pB2[ t, s, xxxx] <- pB[ yrind.surv[xxxx,t,s], t]
        } # xxxx
        mn.phiFY[ t, s] <- mean( phiFY2[ t, s, 1:surv.end[t,s] ] ) 
        mn.phiA[ t, s] <- mean( phiA2[ t, s, 1:surv.end[t,s] ] )
        mn.phiB[ t, s] <- mean( phiB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiFYB[ t, s] <- mean( psiFYB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiAB[ t, s] <- mean( psiAB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiBA[ t, s] <- mean( psiBA2[ t, s, 1:surv.end[t,s] ] )
        mn.pA[ t, s] <- mean( pA2[ t, s, 1:surv.end[t,s] ] )
        mn.pB[ t, s] <- mean( pB2[ t, s, 1:surv.end[t,s] ] )
      }}
    
    for (i in 1:nind){
      for (t in 1:nyr){
        #Survival
        logit(phiFY[i,t]) <- eta[1, site[i,t],t] + eps[1,t] + 
          lmus[1, site[i,t]] + betas[1]*hacked[i]  # first year
        logit(phiA[i,t]) <- eta[2, site[i,t],t] + eps[2,t] + 
          lmus[2, site[i,t]] +  betas[2]*hacked[i] # nonbreeder
        logit(phiB[i,t]) <- eta[3, site[i,t],t] + eps[3,t] + 
          lmus[3, site[i,t]] + betas[3]*hacked[i] # breeder
        #Recruitment
        logit(psiFYB[i,t]) <- eta[4, site[i,t],t] + eps[4,t] + 
          lmus[4, site[i,t]] + betas[4]*hacked[i] # first year to breeder
        logit(psiAB[i,t]) <- eta[5, site[i,t],t] + eps[5,t] + 
          lmus[5, site[i,t]] + betas[5]*hacked[i] # nonbreeder to breeder
        logit(psiBA[i,t]) <- #eta[6, site[i,t],t] + eps[6,t] + 
          lmus[6, site[i,t]] #+ betas[6]*hacked[i] # breeder to nonbreeder
        #Re-encounter
        logit(pA[i,t]) <- eta[7, site[i,t],t] + eps[7,t] + 
          lmus[7, site[i,t]] + betas[7]*hacked[i] # resight of nonbreeders
        logit(pB[i,t]) <- eta[8, site[i,t],t] + eps[8,t] + 
          lmus[8, site[i,t]] + betas[8]*hacked[i] # resight of breeders
      }#t
    }#i
    
    # Define state-transition and observation matrices
    for (i in 1:nind){  
      # Define probabilities of state S(t+1) given S(t)
      for (t in first[i]:(nyr-1)){
        ps[1,i,t,1] <- 0
        ps[1,i,t,2] <- phiFY[i,t] * (1-psiFYB[i,t]) 
        ps[1,i,t,3] <- phiFY[i,t] * psiFYB[i,t]
        ps[1,i,t,4] <- (1-phiFY[i,t])
        
        ps[2,i,t,1] <- 0
        ps[2,i,t,2] <- phiA[i,t] * (1-psiAB[i,t])
        ps[2,i,t,3] <- phiA[i,t] * psiAB[i,t]
        ps[2,i,t,4] <- (1-phiA[i,t])
        
        ps[3,i,t,1] <- 0
        ps[3,i,t,2] <- phiB[i,t] * psiBA[i,t]
        ps[3,i,t,3] <- phiB[i,t] * (1-psiBA[i,t])
        ps[3,i,t,4] <- (1-phiB[i,t])
        
        ps[4,i,t,1] <- 0
        ps[4,i,t,2] <- 0
        ps[4,i,t,3] <- 0
        ps[4,i,t,4] <- 1
        
        # Define probabilities of O(t) given S(t)
        po[1,i,t,1] <- 1 
        po[1,i,t,2] <- 0
        po[1,i,t,3] <- 0
        po[1,i,t,4] <- 0
        
        po[2,i,t,1] <- 0
        po[2,i,t,2] <- pA[i,t]
        po[2,i,t,3] <- 0
        po[2,i,t,4] <- (1-pA[i,t])
        
        po[3,i,t,1] <- 0
        po[3,i,t,2] <- 0
        po[3,i,t,3] <- pB[i,t]
        po[3,i,t,4] <- (1-pB[i,t])
        
        po[4,i,t,1] <- 0
        po[4,i,t,2] <- 0
        po[4,i,t,3] <- 0
        po[4,i,t,4] <- 1
      } #t
    } #i
    
    # Likelihood 
    for (i in 1:nind){
      # Define latent state at first capture
      z[i,first[i]] <- y[i,first[i]]
      for (t in (first[i]+1):nyr){
        # State process: draw S(t) given S(t-1)
        z[i,t] ~ dcat(ps[z[i,t-1], i, t-1, 1:4])
        # Observation process: draw O(t) given S(t)
        y[i,t] ~ dcat(po[z[i,t], i, t-1, 1:4])
      } #t
    } #i
  } )

#**********************
#* Function to run model in NIMBLE
#**********************
run_ipm <- function(info, datl, constl, code){
  library('nimble')
  library('coda')
  # helper function for multivariate normal
  uppertri_mult_diag <- nimbleFunction(
    run = function(mat = double(2), vec = double(1)) {
      returnType(double(2))
      p <- length(vec)
      out <- matrix(nrow = p, ncol = p, init = FALSE)
      for(i in 1:p){
        out[ ,i] <- mat[ ,i] * vec[i]
      }
      return(out)
    })
  assign('uppertri_mult_diag', uppertri_mult_diag, envir = .GlobalEnv)
  
  params <- c(# pop growth 
    #"lambda",
    # fecundity
    "lmu.f", "gamma", "rr", "mn.f", 
    # survival 
    "mus", "lmus", "betas",
    # abundance
    # "NB", "NF", "NFY", "N",
    # "r",
    # "N", "Ntot",
    # error terms
    "eps", "sds", "Ustar", "U", "R",
    "eta", "sds2", "Ustar2", "U2", "R2",
    # yearly summaries
    'mn.phiFY', 'mn.phiA', 'mn.phiB',
    'mn.psiFYB', 'mn.psiAB', 'mn.psiBA',
    'mn.pA', 'mn.pB',
    # goodness of fit
    "f.dmape.obs", "f.dmape.rep",
    "f.tvm.obs", "f.tvm.rep"
    # "dmape.obs", "dmape.rep",
    # "tvm.obs", "tvm.rep",
    # "tturn.obs", "tturn.rep"
  )
  
  #n.chains=1; n.thin=200; n.iter=500000; n.burnin=300000
  #n.chains=1; n.thin=50; n.iter=100000; n.burnin=50000
  #n.chains=1; n.thin=10; n.iter=20000; n.burnin=10000
  n.chains=1; n.thin=200; n.iter=500000; n.burnin=300000
  
  mod <- nimbleModel(code, 
                     constants = constl, 
                     data = datl, 
                     inits = info$inits, 
                     buildDerivs = FALSE, # doesn't work when TRUE, no hope for HMC
                     calculate=T 
  ) 
  
  cmod <- compileNimble(mod, showCompilerOutput = TRUE)
  confhmc <- configureMCMC(mod)
  confhmc$setMonitors(params)
  hmc <- buildMCMC(confhmc)
  chmc <- compileNimble(hmc, project = mod, 
                        resetFunctions = TRUE,
                        showCompilerOutput = TRUE)
  
  post <- runMCMC(chmc,
                  niter = n.iter, 
                  nburnin = n.burnin,
                  nchains = n.chains,
                  thin = n.thin,
                  samplesAsCodaMCMC = T, 
                  setSeed = info$seed)
  
  return(post)
} # run_ipm function end

#*****************
#* Run chains in parallel
#*****************
this_cluster <- makeCluster(10)
post <- parLapply(cl = this_cluster, 
                  X = par_info, 
                  fun = run_ipm, 
                  datl = datl, 
                  constl = constl, 
                  code = mycode)
stopCluster(this_cluster)

save(post, mycode,
     file="/bsuscratch/brianrolek/riha_ipm/outputs/ipm_simp.rdata")

# save(post, mycode,
#      file="C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm_statespace.rdata")

1.2 Integrated Population Model

library('nimble')
library('parallel')
library ('coda')
load("/bsuscratch/brianrolek/riha_ipm/data.rdata")

#**********************
#* Parameter descriptions
#**********************
# Mark-resight-recovery data
#   Observations (po) = y  
#     1 seen first-year (age=0, just before 1st b-day)
#     2 seen nonbreeder
#     3 seen breeder
#     4 not seen
#   States (ps)
#     1 alive first-year
#     2 alive nonbreeder
#     3 alive breeder
#     4 dead
#   Groups
#     1 wild-born
#     2 translocated and hacked
###################################################
# PARAMETERS
#   phiFY: survival probability first year 
#   phiA: survival probability nonbreeders
#   phiB: survival probability breeders
#   psiFYB: recruitment probability from first-year to breeder
#   psiAB: recruitment probability from nonbreeders to breeder
#   psiBA: recruitment probability from breeder to nonbreeders 
#   pFY: resight probability first-years
#   pA: resight probability nonbreeders
#   pB: resight probability breeders
#   lmu.prod: mean productivity (males and females) per territory on the log scale
#   sds: standard deviations for multivariate normal random effects for time and site
#   R: correlation coefficients for multivariate normal random effects for time and site
#   lambda: population growth rate (derived)
#   extinct: binary indicator of extirpation at a site
#   gamma: coefficient of effect from nest treatments
#   betas: coefficient of effect from translocations
#   deltas: coefficient of effect from survey effort
#   mus: overall means for survival, recruitment, and detections
#   r and rr: "r" parameter for negative binomial distribution
#             also described as omega in manuscript

#**********************
#* Model code
#**********************
mycode <- nimbleCode(
  {
    ###################################################
    # Priors and constraints
    ###################################################
    # survival, recruitment, and detection can be correlated
    for (k in 1:8){
      betas[k] ~ dunif(-20, 20)  # prior for translocations coefficients
      deltas[k] ~ dunif(-20, 20) # prior for survey effort coefficients
    } # k
    
    for (j in 1:8){
      for (s in 1:nsite){
        lmus[j,s] <- logit(mus[j,s])
        mus[j,s] ~ dbeta(1,1) # prior for overall means
      }}     # 
    
    # Temporal random effects and correlations among all sites
    for (j in 1:p){ sds[j] ~ dexp(1) }# prior for temporal variation
    # estimated using the multivariate normal distribution
    R[1:p,1:p] <- t(Ustar[1:p,1:p]) %*% Ustar[1:p,1:p] # calculate rhos, correlation coefficients
    Ustar[1:p,1:p] ~ dlkj_corr_cholesky(eta=1.1, p=p) # Ustar is the Cholesky decomposition of the correlation matrix
    U[1:p,1:p] <- uppertri_mult_diag(Ustar[1:p, 1:p], sds[1:p])
    # multivariate normal for temporal variance
    for (t in 1:nyr){ # survival params only have nyr-1, no problem to simulate from however
      eps[1:p,t] ~ dmnorm(mu.zeroes[1:p],
                          cholesky = U[1:p, 1:p], prec_param = 0)
    }
    
    # Temporal random effects and correlations between sites
    for (jj in 1:p2){ sds2[jj] ~ dexp(1) }# prior for temporal variation
    # estimated using the multivariate normal distribution
    R2[1:p2,1:p2] <- t(Ustar2[1:p2,1:p2]) %*% Ustar2[1:p2,1:p2] # calculate rhos, correlation coefficients
    Ustar2[1:p2,1:p2] ~ dlkj_corr_cholesky(eta=1.1, p=p2) # Ustar is the Cholesky decomposition of the correlation matrix
    U2[1:p2,1:p2] <- uppertri_mult_diag(Ustar2[1:p2, 1:p2], sds2[1:p2])
    # multivariate normal for temporal variance
    for (t in 1:nyr){ # survival params only have nyr-1, no problem to simulate from however
      for (s in 1:nsite){
        eta[1:p2,s,t] ~ dmnorm(mu.zeroes2[1:p2],
                               cholesky = U2[1:p2, 1:p2], prec_param = 0)
      } } # s t 
    #######################
    # Derived params
    #######################
    for (s in 1:nsite){
      for (t in 1:(nyr-1)){
        lambda[t, s] <-  Ntot[t+1, s]/(Ntot[t, s])
        loglambda[t, s] <- log(lambda[t, s])
      }} #t
    
    ###############################
    # Likelihood for productivity
    ###############################
    # Priors
    for (s in 1:nsite){
      lmu.prod[s] ~ dnorm(0, sd=5)
    } # s
    gamma ~ dunif(-20, 20)
    rr ~ dexp(0.05)
    
    # Productivity likelihood      
    for (k in 1:npairsobs){
      prod[k] ~ dnegbin(ppp[k], rr)
      ppp[k] <- rr/(rr+mu.prod[k])
      log(mu.prod[k]) <- lmu.prod[site.pair[k]] +  
        gamma*treat.pair[k] + 
        eps[9, year.pair[k] ] + 
        eta[9, site.pair[k], year.pair[k] ] 
    } # k
    # Derive yearly productivity for population model
    # need to reorder because nimble doesn't 
    # handle nonconsecutive indices
    # yrind.pair is a matrix of indices for each site
    for (t in 1:nyr){
      for (s in 1:nsite){
        for (xxx in 1:pair.end[t,s]){
          prodmat[t,s,xxx] <- mu.prod[ yrind.pair[xxx,t,s] ]
        } # xxx
        mn.prod[t,s] <- mean( prodmat[t,s,1:pair.end[t,s]] )
      }} # s t
    
    # GOF for productivity
    for (k in 1:npairsobs){
      f.obs[k] <- prod[k] # observed counts
      f.exp[k] <- mu.prod[k] # expected counts adult breeder
      f.rep[k] ~ dnegbin(ppp[k], rr) # expected counts
      f.dssm.obs[k] <- abs( ( f.obs[k] - f.exp[k] ) / (f.obs[k]+0.001) )
      f.dssm.rep[k] <- abs( ( f.rep[k] - f.exp[k] ) / (f.rep[k]+0.001) )
    } # k
    f.dmape.obs <- sum(f.dssm.obs[1:npairsobs])
    f.dmape.rep <- sum(f.dssm.rep[1:npairsobs])
    f.tvm.obs <- sd(brood[1:npairsobs])^2/mean(brood[1:npairsobs])
    f.tvm.rep <- sd(f.rep[1:npairsobs])^2/mean(f.rep[1:npairsobs])
    
    ################################
    # Likelihood for counts
    ################################
    # Abundance for year=1
    for (v in 1:7){ 
      for (s in 1:nsite){
        # subtract one to allow dcat to include zero
        N[v, 1, s] <- N2[v, 1, s] - 1 
        N2[v, 1, s] ~ dcat(pPrior[v, 1:s.end[v,s], s]) # Priors differ for FYs and adults
      }} # s t
    # Abundance for years > 1
    for (t in 1:(nyr-1)){
      for (s in 1:nsite){
        # Number of wild born juvs
        N[1, t+1, s] ~ dpois( (NFY[t, s]*mn.phiFY[t, s]*mn.psiFYB[t, s] + # first year breeders
                                 NF[t, s]*mn.phiA[t, s]*mn.psiAB[t, s] + # nonbreeders to breeders
                                 NB[t, s]*mn.phiB[t, s]*(1-mn.psiBA[t, s])) # breeders remaining
                              *mn.prod[t+1, s]/2 ) # end Poisson
        # Abundance of nonbreeders
        N[2, t+1, s] ~ dbin(mn.phiFY[t, s]*(1-mn.psiFYB[t, s]), NFY[t, s]) # Nestlings to nonbreeders
        N[3, t+1, s] ~ dbin(mn.phiA[t, s]*(1-mn.psiAB[t, s]), NF[t, s]) # Nonbreeders to nonbreeders
        N[4, t+1, s] ~ dbin(mn.phiB[t, s]*mn.psiBA[t, s], NB[t, s]) # Breeders to nonbreeders
        # Abundance of breeders
        N[5, t+1, s] ~ dbin(mn.phiFY[t, s]*mn.psiFYB[t, s], NFY[t, s]) # Nestlings to second year breeders
        N[6, t+1, s] ~ dbin(mn.phiA[t, s]*mn.psiAB[t, s], NF[t, s]) # Nonbreeder to breeder
        N[7, t+1, s] ~ dbin(mn.phiB[t, s]*(1-mn.psiBA[t, s]), NB[t, s]) # Breeder to breeder
      }} # s t
    
    # Count likelihoods, state-space model, and observation process    
    for (t in 1:nyr){
      for (s in 1:nsite){
        countsAdults[t, s] ~ dpois(lamAD[t, s]) # adult females 
        constraint_data[t, s] ~ dconstraint( (N[1, t, s] + hacked.counts[t, s]) >= 0 ) # constrain N1+hacked.counts to be >=0
        NFY[t, s] <- N[1, t, s] + hacked.counts[t, s] # Transfers translocated first-year females
        NF[t, s] <- sum(N[2:4, t, s]) # number of adult nonbreeder females
        NB[t, s] <- sum(N[5:7, t, s]) # number of adult breeder females
        NAD[t, s] <- NF[t, s] + NB[t, s] # number of adults
        Ntot[t, s] <- sum(N[1:7, t, s]) # total number of females
        log(lamAD[t, s]) <- log(NAD[t, s]) + deltas[1]*effort2[t, s] + deltas[2]*effort2[t, s]^2
        log(lamFY[t, s]) <- log(N[1, t, s]) + deltas[3]*effort2[t, s] + deltas[4]*effort2[t, s]^2
        }# s
      # First-years at different sites have different distributions
      # for better model fit
      countsFY[t, 1] ~ dpois(lamFY[t, 1]) # doesn't have any zeroes so poisson
      countsFY[t, 2] ~ dnegbin(pp[t], r) # first year females, includes translocated/hacked
      pp[t] <- r/(r+(lamFY[t, 2] ))
    } # t
    r ~ dexp(0.05)
    ###################
    # Assess GOF of the state-space models for counts
    # Step 1: Compute statistic for observed data
    # Step 2: Use discrepancy measure: mean absolute error
    # Step 3: Use test statistic: number of turns
    ###################
    for (t in 1:nyr){
      c.repFY[t, 1] ~ dpois( lamFY[t, 1] )
      c.repFY[t, 2] ~ dnegbin( pp[t], r )
      for (s in 1:nsite){
        c.expAD[t, s] <- lamAD[t, s]  # expected counts adult breeder
        c.expFY[t, s] <- lamFY[t, s]
        c.obsAD[t, s] <- countsAdults[t, s]
        c.obsFY[t, s] <- countsFY[t, s]  # first year
        c.repAD[t, s] ~ dpois( lamAD[t, s] ) # simulated counts
        dssm.obsAD[t, s] <- abs( ( (c.obsAD[t, s]) - (c.expAD[t, s]) ) / (c.obsAD[t, s]+0.001)  )
        dssm.obsFY[t, s] <- abs( ( (c.obsFY[t, s]) - (c.expFY[t, s]) ) / (c.obsFY[t, s]+0.001)  )
        dssm.repAD[t, s] <- abs( ( (c.repAD[t, s]) - (c.expAD[t, s]) ) / (c.repAD[t, s]+0.001) )
        dssm.repFY[t, s] <- abs( ( (c.repFY[t, s]) - (c.expFY[t, s]) ) / (c.repFY[t, s]+0.001) )
      }} # t
    dmape.obs[1] <- sum(dssm.obsAD[1:nyr, 1:nsite])
    dmape.obs[2] <- sum(dssm.obsFY[1:nyr, 1:nsite])
    dmape.rep[1] <- sum(dssm.repAD[1:nyr, 1:nsite])
    dmape.rep[2] <- sum(dssm.repFY[1:nyr, 1:nsite])
    tvm.obs[1] <- sd(c.obsB[1:nyr, 1:nsite])^2/mean(c.obsB[1:nyr, 1:nsite])
    tvm.obs[2] <- sd(c.obsFY[1:nyr, 1:nsite])^2/mean(c.obsFY[1:nyr, 1:nsite])
    tvm.rep[1] <- sd(c.repB[1:nyr, 1:nsite])^2/mean(c.repB[1:nyr, 1:nsite])
    tvm.rep[2] <- sd(c.repFY[1:nyr, 1:nsite])^2/mean(c.repFY[1:nyr, 1:nsite])
    # # Test statistic for number of turns
    for (s in 1:nsite){
      for (t in 1:(nyr-2)){
        tt1.obsAD[t, s] <- step(c.obsAD[t+2, s] - c.obsAD[t+1, s])
        tt2.obsAD[t, s] <- step(c.obsAD[t+1, s] - c.obsAD[t, s])
        tt3.obsAD[t, s] <- equals(tt1.obsAD[t, s] + tt2.obsAD[t, s], 1)
        tt1.obsFY[t, s] <- step(c.obsFY[t+2, s] - c.obsFY[t+1, s])
        tt2.obsFY[t, s] <- step(c.obsFY[t+1, s] - c.obsFY[t, s])
        tt3.obsFY[t, s] <- equals(tt1.obsFY[t, s] + tt2.obsFY[t, s], 1)
      }} # t
    tturn.obs[1] <- sum(tt3.obsAD[1:(nyr-2), 1:nsite])
    tturn.obs[2] <- sum(tt3.obsFY[1:(nyr-2), 1:nsite])
    for (s in 1:nsite){
      for (t in 1:(nyr-2)){
        tt1.repAD[t, s] <- step(c.repAD[t+2, s] - c.repAD[t+1, s])
        tt2.repAD[t, s] <- step(c.repAD[t+1, s] - c.repAD[t, s])
        tt3.repAD[t, s] <- equals(tt1.repAD[t, s] + tt2.repAD[t, s], 1)
        tt1.repFY[t, s] <- step(c.repFY[t+2, s] - c.repFY[t+1, s])
        tt2.repFY[t, s] <- step(c.repFY[t+1, s] - c.repFY[t, s])
        tt3.repFY[t, s] <- equals(tt1.repFY[t, s] + tt2.repFY[t, s], 1)
      }} # t
    tturn.rep[1] <- sum(tt3.repAD[1:(nyr-2), 1:nsite])
    tturn.rep[2] <- sum(tt3.repFY[1:(nyr-2), 1:nsite])
    
    ################################
    # Likelihood for survival
    ################################ 
    # Calculate averages for sites each year for integration
    for (t in 1:nyr){
      for (s in 1:nsite){
        for (xxxx in 1:surv.end[t,s]){
          # Reorder because nimble doesn't 
          # handle nonconsecutive indices
          # yrind.surv is a matrix of indices for each site
          phiFY2[ t, s, xxxx] <- phiFY[ yrind.surv[xxxx,t,s], t]
          phiA2[ t, s, xxxx] <- phiA[ yrind.surv[xxxx,t,s], t]
          phiB2[ t, s, xxxx] <- phiB[ yrind.surv[xxxx,t,s], t]
          psiFYB2[ t, s, xxxx] <- psiFYB[ yrind.surv[xxxx,t,s], t]
          psiAB2[ t, s, xxxx] <- psiAB[ yrind.surv[xxxx,t,s], t]
          psiBA2[ t, s, xxxx] <- psiBA[ yrind.surv[xxxx,t,s], t]
          pA2[ t, s, xxxx] <- pA[ yrind.surv[xxxx,t,s], t]
          pB2[ t, s, xxxx] <- pB[ yrind.surv[xxxx,t,s], t]
        } # xxxx
        mn.phiFY[ t, s] <- mean( phiFY2[ t, s, 1:surv.end[t,s] ] ) 
        mn.phiA[ t, s] <- mean( phiA2[ t, s, 1:surv.end[t,s] ] )
        mn.phiB[ t, s] <- mean( phiB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiFYB[ t, s] <- mean( psiFYB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiAB[ t, s] <- mean( psiAB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiBA[ t, s] <- mean( psiBA2[ t, s, 1:surv.end[t,s] ] )
        mn.pA[ t, s] <- mean( pA2[ t, s, 1:surv.end[t,s] ] )
        mn.pB[ t, s] <- mean( pB2[ t, s, 1:surv.end[t,s] ] )
      }}
    
    for (i in 1:nind){
      for (t in 1:nyr){
        #Survival
        logit(phiFY[i,t]) <- eta[1, site[i,t],t] + eps[1,t] + 
          lmus[1, site[i,t]] + betas[1]*hacked[i]  # first year
        logit(phiA[i,t]) <- eta[2, site[i,t],t] + eps[2,t] + 
          lmus[2, site[i,t]] +  betas[2]*hacked[i] # nonbreeder
        logit(phiB[i,t]) <- eta[3, site[i,t],t] + eps[3,t] + 
          lmus[3, site[i,t]] + betas[3]*hacked[i] # breeder
        #Recruitment
        logit(psiFYB[i,t]) <- eta[4, site[i,t],t] + eps[4,t] + 
          lmus[4, site[i,t]] + betas[4]*hacked[i] # first year to breeder
        logit(psiAB[i,t]) <- eta[5, site[i,t],t] + eps[5,t] + 
          lmus[5, site[i,t]] + betas[5]*hacked[i] # nonbreeder to breeder
        logit(psiBA[i,t]) <- #eta[6, site[i,t],t] + eps[6,t] + 
          lmus[6, site[i,t]] #+ betas[6]*hacked[i] # breeder to nonbreeder
        #Re-encounter
        logit(pA[i,t]) <- eta[7, site[i,t],t] + eps[7,t] + 
          lmus[7, site[i,t]] + betas[7]*hacked[i] +
          deltas[5]*effort2[t, site[i,t]] + deltas[6]*effort2[t, site[i,t]]^2# resight of nonbreeders
        logit(pB[i,t]) <- eta[8, site[i,t],t] + eps[8,t] + 
          lmus[8, site[i,t]] + betas[8]*hacked[i] + 
          deltas[7]*effort2[t, site[i,t]] + deltas[8]*effort2[t, site[i,t]]^2# resight of breeders
      }#t
    }#i
    
    # Define state-transition and observation matrices
    for (i in 1:nind){  
      # Define probabilities of state S(t+1) given S(t)
      for (t in first[i]:(nyr-1)){
        ps[1,i,t,1] <- 0
        ps[1,i,t,2] <- phiFY[i,t] * (1-psiFYB[i,t]) 
        ps[1,i,t,3] <- phiFY[i,t] * psiFYB[i,t]
        ps[1,i,t,4] <- (1-phiFY[i,t])
        
        ps[2,i,t,1] <- 0
        ps[2,i,t,2] <- phiA[i,t] * (1-psiAB[i,t])
        ps[2,i,t,3] <- phiA[i,t] * psiAB[i,t]
        ps[2,i,t,4] <- (1-phiA[i,t])
        
        ps[3,i,t,1] <- 0
        ps[3,i,t,2] <- phiB[i,t] * psiBA[i,t]
        ps[3,i,t,3] <- phiB[i,t] * (1-psiBA[i,t])
        ps[3,i,t,4] <- (1-phiB[i,t])
        
        ps[4,i,t,1] <- 0
        ps[4,i,t,2] <- 0
        ps[4,i,t,3] <- 0
        ps[4,i,t,4] <- 1
        
        # Define probabilities of O(t) given S(t)
        po[1,i,t,1] <- 1 
        po[1,i,t,2] <- 0
        po[1,i,t,3] <- 0
        po[1,i,t,4] <- 0
        
        po[2,i,t,1] <- 0
        po[2,i,t,2] <- pA[i,t]
        po[2,i,t,3] <- 0
        po[2,i,t,4] <- (1-pA[i,t])
        
        po[3,i,t,1] <- 0
        po[3,i,t,2] <- 0
        po[3,i,t,3] <- pB[i,t]
        po[3,i,t,4] <- (1-pB[i,t])
        
        po[4,i,t,1] <- 0
        po[4,i,t,2] <- 0
        po[4,i,t,3] <- 0
        po[4,i,t,4] <- 1
      } #t
    } #i
    
    # Likelihood 
    for (i in 1:nind){
      # Define latent state at first capture
      z[i,first[i]] <- y[i,first[i]]
      for (t in (first[i]+1):nyr){
        # State process: draw S(t) given S(t-1)
        z[i,t] ~ dcat(ps[z[i,t-1], i, t-1, 1:4])
        # Observation process: draw O(t) given S(t)
        y[i,t] ~ dcat(po[z[i,t], i, t-1, 1:4])
      } #t
    } #i
  } )

#**********************
#* Function to run model in NIMBLE
#**********************
run_ipm <- function(info, datl, constl, code){
  library('nimble')
  library('coda')
  # helper function for multivariate normal
  uppertri_mult_diag <- nimbleFunction(
    run = function(mat = double(2), vec = double(1)) {
      returnType(double(2))
      p <- length(vec)
      out <- matrix(nrow = p, ncol = p, init = FALSE)
      for(i in 1:p){
        out[ ,i] <- mat[ ,i] * vec[i]
      }
      return(out)
    })
  assign('uppertri_mult_diag', uppertri_mult_diag, envir = .GlobalEnv)
  
  params <- c(
    # pop growth 
    "lambda",
    # fecundity
    "lmu.prod", "gamma", "rr", "mn.prod", 
    # survival 
    "mus", "lmus", "betas", "deltas",
    # abundance
    "NB", "NF", "NFY", "N", "NAD",
    "r",
    "N", "Ntot",
    # error terms
    "eps", "sds", "Ustar", "U", "R",
    "eta", "sds2", "Ustar2", "U2", "R2",
    # yearly summaries
    'mn.phiFY', 'mn.phiA', 'mn.phiB',
    'mn.psiFYB', 'mn.psiAB', 'mn.psiBA',
    'mn.pA', 'mn.pB',
    # goodness of fit
    "f.dmape.obs", "f.dmape.rep",
    "f.tvm.obs", "f.tvm.rep",
    "dmape.obs", "dmape.rep",
    "tvm.obs", "tvm.rep",
    "tturn.obs", "tturn.rep"
  )
  n.chains=1; n.thin=1; n.iter=500; n.burnin=100
  #n.chains=1; n.thin=200; n.iter=600000; n.burnin=400000
  
  mod <- nimbleModel(code, 
                     constants = constl, 
                     data = datl, 
                     inits = info$inits, 
                     buildDerivs = FALSE, # doesn't work when TRUE, no hope for HMC
                     calculate=T 
  ) 
  
  cmod <- compileNimble(mod, showCompilerOutput = TRUE)
  confhmc <- configureMCMC(mod)
  confhmc$setMonitors(params)
  hmc <- buildMCMC(confhmc)
  chmc <- compileNimble(hmc, project = mod, 
                        resetFunctions = TRUE,
                        showCompilerOutput = TRUE)
  
  post <- runMCMC(chmc,
                  niter = n.iter, 
                  nburnin = n.burnin,
                  nchains = n.chains,
                  thin = n.thin,
                  samplesAsCodaMCMC = T, 
                  setSeed = info$seed)
  
  return(post)
} # run_ipm function end

#*****************
#* Run chains in parallel
#*****************
this_cluster <- makeCluster(20)
post <- parLapply(cl = this_cluster, 
                  X = par_info, 
                  fun = run_ipm, 
                  datl = datl, 
                  constl = constl, 
                  code = mycode)
stopCluster(this_cluster)
# save(post, mycode,
#      file="/bsuscratch/brianrolek/riha_ipm/outputs/ipm_longrun.rdata")

1.3 Population Viability Analysis

library('nimble')
library('parallel')
library ('coda')
load("/bsuscratch/brianrolek/riha_ipm/data.rdata")
source("/bsuscratch/brianrolek/riha_ipm/MCMCvis.R")
constl$K <- 100
cpus <- 10
#**********************
#* Extend "data" needed for PVA projections and scenarios
#**********************
constl$SC <- 12

hacked.counts <- array(0, dim=c(constl$SC, constl$nyr+constl$K, 3))
for (sc in 1:constl$SC){
  hacked.counts[sc,1:constl$nyr,1:3] <- constl$hacked.counts 
  hacked.counts[sc,14:(constl$nyr+constl$K),1:3] <- 
    c(0, 0, 0, 1, 1, 1)[sc]*cbind(rep(-10, constl$K), rep(10, constl$K), rep(0, constl$K) )
}
constl$hacked.counts <- hacked.counts
datl$constraint_data <- rbind(datl$constraint_data, array(1, dim=c(constl$K,2)) )
#constl$treat.pair2 <- c(1, 0, 1, 1, 0, 1) # treated or not
constl$treat.pair2 <- c(0, 1, 1, 0, 1, 1) # treated or not
constl$hacked2 <- c(0, 0, 0, 1, 1, 1) # scenario- hacked or not
constl$effort2 <- rbind(constl$effort2, array(0, dim=c(constl$K,2), dimnames=list(2024:(2024+constl$K-1), c("LH", "PC"))))
# set number of treated nests for all scenarios
constl$num.treated <- c(0, 10, 100, 0, 10, 100, 0, 10, 100, 0, 10, 100) # 100s sub in for All but are over-ridden in model code
constl$treat.all <- c(0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1)
#**********************
#* Specify priors taken from posterior of IPM
#**********************
#* This helps ensure that all chains run
load("/bsuscratch/brianrolek/riha_ipm/outputs/ipm_longrun.rdata")
# Identify chains with NAs that failed to initialize
out <- lapply(post, as.mcmc)
NAlist <- c()
for (i in 1:length(out)){
  NAlist[i] <- any (is.na(out[[i]][,1:286]) | out[[i]][,1:286]<0)
}
out <- out[!NAlist]
outp <- MCMCpstr(out, type="chains")
!NAlist

mn.lam <- apply(outp$lambda, c(2,3), mean)
lam <- apply(mn.lam, 1, mean)

Ni.func <- function (){
  Ni <- outp$N[1:7,1:constl$nyr,1:2,
               sample(seq(1, 5200, by=400), 1, replace = F)]
  Ni.pva <- array(NA, dim=c(constl$SC, dim(Ni)+c(0,constl$K,0)))
  for (sc in 1:constl$SC){
    Ni.pva[sc, 1:7, 1:constl$nyr, 1:2] <- Ni
    for (t in constl$nyr:(constl$nyr+constl$K)){
      for (s in 1:2){
      Ni.pva[sc, 1:7, t, s] <- Ni[1:7, 13, s] #round(Ni[1:7, 13, s]*lam[s]^(t-13))
    }}} # t sc
  return(Ni.pva)
} # function

u1 <- apply(outp$Ustar, c(1,2), mean)
u2 <- apply(outp$Ustar2, c(1,2), mean)

inits.func.pva <- function (){
  list(  
    # fecundity 
    lmu.prod = apply(outp$lmu.prod, 1, mean),
    gamma = mean(outp$gamma), 
    rr = mean(outp$rr),
    # survival
    z = par_info[[1]]$inits$z, 
    mus = apply(outp$mus, c(1,2), mean), 
    betas = apply(outp$betas, 1, mean),
    deltas = apply(outp$deltas, 1, mean),
    sds = apply(outp$sds, 1, mean),
    Ustar = u1,
    sds2 =  apply(outp$sds2, 1, mean),
    Ustar2 = u1,
    # counts
    countsAdults= matrix(c(374, 335, 305, 295, rep(NA, length(2015:2023)), rep(NA, length(2011:2023)) ), nrow=13), 
    r = mean(outp$r),
    N = Ni.func()
  )}

# Set seed then save for reproducibility
# then randomly draw for each chain
set.seed(1)
seeds <- sample(1:1000000, size=cpus, replace=FALSE)
par_info_pva <- list()
for (i in 1:cpus){
  par_info_pva[[i]] <- list(seed=seeds[i], inits = inits.func.pva())
}

#**********************
#* Parameter descriptions
#**********************
# Mark-resight-recovery data
#   Observations (po) = y  
#     1 seen first-year (age=0, just before 1st b-day)
#     2 seen nonbreeder
#     3 seen breeder
#     4 not seen
#   States (ps)
#     1 alive first-year
#     2 alive nonbreeder
#     3 alive breeder
#     4 dead
#   Groups
#     1 wild-born
#     2 translocated and hacked
###################################################
# PARAMETERS
#   phiFY: survival probability first year 
#   phiA: survival probability nonbreeders
#   phiB: survival probability breeders
#   psiFYB: recruitment probability from first-year to breeder
#   psiAB: recruitment probability from nonbreeders to breeder
#   psiBA: recruitment probability from breeder to nonbreeders 
#   pFY: resight probability first-years
#   pA: resight probability nonbreeders
#   pB: resight probability breeders
#   lmu.prod: mean productivity (males and females) per territory on the log scale
#   sds: standard deviations for multivariate normal random effects for time and site
#   R: correlation coefficients for multivariate normal random effects for time and site
#   lambda: population growth rate (derived)
#   extinct: binary indicator of extirpation at a site
#   gamma: coefficient of effect from nest treatments
#   betas: coefficient of effect from translocations
#   deltas: coefficient of effect from survey effort
#   mus: overall means for survival, recruitment, and detections
#   r and rr: "r" parameter for negative binomial distribution
#             also described as omega in manuscript

#**********************
#* Model code
#**********************
mycode <- nimbleCode(
  {
    ###################################################
    # Priors and constraints
    ###################################################
    # survival, recruitment, and detection can be correlated
    for (k in 1:8){ 
      betas[k] ~ dunif(-20, 20)  # prior for coefficients
      deltas[k] ~ dunif(-20, 20)
    } # k
    
    for (j in 1:8){
      for (s in 1:nsite){
        lmus[j,s] <- logit(mus[j,s])
        mus[j,s] ~ dbeta(1,1) # prior for means
      }}     # m population #s sex #h hacked
    
    # Temporal random effects and correlations among all sites
    for (j in 1:p){ sds[j] ~ dexp(1) }# prior for temporal variation
    # estimated using the multivariate normal distribution
    R[1:p,1:p] <- t(Ustar[1:p,1:p]) %*% Ustar[1:p,1:p] # calculate rhos, correlation coefficients
    Ustar[1:p,1:p] ~ dlkj_corr_cholesky(eta=1.1, p=p) # Ustar is the Cholesky decomposition of the correlation matrix
    U[1:p,1:p] <- uppertri_mult_diag(Ustar[1:p, 1:p], sds[1:p])
    # multivariate normal for temporal variance
    for (t in 1:(nyr+K)){ # survival params only have nyr-1, no problem to simulate from however
      eps[1:p,t] ~ dmnorm(mu.zeroes[1:p],
                          cholesky = U[1:p, 1:p], prec_param = 0)
    }
    
    # Temporal random effects and correlations between sites
    for (jj in 1:p2){ sds2[jj] ~ dexp(1) }# prior for temporal variation
    # estimated using the multivariate normal distribution
    R2[1:p2,1:p2] <- t(Ustar2[1:p2,1:p2]) %*% Ustar2[1:p2,1:p2] # calculate rhos, correlation coefficients
    Ustar2[1:p2,1:p2] ~ dlkj_corr_cholesky(eta=1.1, p=p2) # Ustar is the Cholesky decomposition of the correlation matrix
    U2[1:p2,1:p2] <- uppertri_mult_diag(Ustar2[1:p2, 1:p2], sds2[1:p2])
    # multivariate normal for temporal variance
    for (t in 1:(nyr+K) ){ # survival params only have nyr-1, no problem to simulate from however
      for (s in 1:nsite){
        eta[1:p2,s,t] ~ dmnorm(mu.zeroes2[1:p2],
                               cholesky = U2[1:p2, 1:p2], prec_param = 0)
      } } # s t 
    #######################
    # Derived params
    #######################
    for(sc in 1:SC){
    for (s in 1:nsite){
      for (t in 1:(nyr-1+K)){
        lambda[sc, t, s] <-  Ntot[sc, t+1, s]/(Ntot[sc, t, s])
        loglambda[sc, t, s] <- log(lambda[sc, t, s])
      }} #t
    
    for (s in 1:nsite){
      for (t in 1:K){
        extinct[sc, t, s] <- equals(Ntot[sc, nyr+t, s], 0)
        extinctAD[sc, t, s] <- equals(NAD[sc, nyr+t, s], 0)
        extinctB[sc, t, s] <- equals(NB[sc, nyr+t, s], 0)
      }}
    } # sc
    
    ###############################
    # Likelihood for productivity
    ###############################
    # Priors 
    for (s in 1:nsite){
      lmu.prod[s] ~ dnorm(0, sd=5)
    } # s
    gamma ~ dunif(-20, 20)
    rr ~ dexp(0.05)
    
    # Productivity likelihood       
    for (k in 1:npairsobs){
      prod[k] ~ dnegbin(ppp[k], rr)
      ppp[k] <- rr/(rr+mu.prod[k])
      log(mu.prod[k]) <- lmu.prod[site.pair[k]] +  
        gamma*treat.pair[k] + 
        eps[9, year.pair[k] ] + 
        eta[9, site.pair[k], year.pair[k] ] 
    } # k
    # Derive yearly productivity for population model
    # need to reorder because nimble doesn't 
    # handle nonconsecutive indices
    # yrind.pair is a matrix of indices for each site
    for (t in 1:nyr){
      for (s in 1:nsite){
        for (xxx in 1:pair.end[t,s]){
          prodmat[t,s,xxx] <- mu.prod[ yrind.pair[xxx,t,s] ]
        } # xxx
        mn.prod[1,t,s] <- mean( prodmat[t,s,1:pair.end[t,s]] )
        for(sc in 2:SC){
          mn.prod[sc, t, s] <- mn.prod[1,t,s]
      }}} # s t
    # Future projections- fecundity
      for(sc in 1:SC){
        for (t in (nyr+1):(nyr+K) ){
          for (s in 1:nsite){
        # Calculate percent of nests treated 
        perc.treat[sc, t, s] <-  treat.all[sc] + # If scenario =3 or 6, prop =1.0
                                step( NB[sc, t, s]-num.treated[sc] ) * (num.treated[sc] + 0.001)/(NB[sc, t, s]+0.001) * (1-treat.all[sc]) + # NB>=10, calculate proportion
                                step( NB[sc, t, s]-1 ) * 1-step( NB[sc, t, s]-num.treated[sc] ) * (1-treat.all[sc])  # NB>1 and NB<10 100%
        log(mn.prod[sc, t, s]) <- lmu.prod[s] +  
                              gamma*treat.pair2[sc]*perc.treat[sc, t, s] + 
                              eps[9, t] + 
                              eta[9, s, t] 
  }} } # s t sc
    
    ################################
    # Likelihood for counts
    ################################
    # Abundance for year=1
      for (v in 1:7){ 
      for (s in 1:nsite){
        # Subtract one to allow dcat to include zero
        N[1, v, 1, s] <- N2[1, v, 1, s] - 1 
        N2[1, v, 1, s] ~ dcat(pPrior[v, 1:s.end[v,s], s])
        # Assign estimates for years with data to all scenarios
        for (sc in 2:SC){
          N[sc, v, 1, s] <- N[1, v, 1, s]
        }
      }} # s t

    # Abundance for years > 1
    for (t in 1:(nyr-1)){
      for (s in 1:nsite){
        # Number of wild born juvs
        N[1, 1, t+1, s] ~ dpois( (NFY[1, t, s]*mn.phiFY[1, t, s]*mn.psiFYB[1, t, s] + # first year breeders
                                 NF[1, t, s]*mn.phiA[1, t, s]*mn.psiAB[1, t, s] + # nonbreeders to breeders
                                 NB[1, t, s]*mn.phiB[1, t, s]*(1-mn.psiBA[1, t, s])) # breeders remaining
                              *(mn.prod[1, t+1, s]/2) ) # end Poisson
        # Abundance of nonbreeders
        N[1, 2, t+1, s] ~ dbin(mn.phiFY[1, t, s]*(1-mn.psiFYB[1, t, s]), NFY[1, t, s]) # Nestlings to second year nonbreeders
        N[1, 3, t+1, s] ~ dbin(mn.phiA[1, t, s]*(1-mn.psiAB[1, t, s]), NF[1, t, s]) # Nonbreeders to nonbreeders
        N[1, 4, t+1, s] ~ dbin(mn.phiB[1, t, s]*mn.psiBA[1, t, s], NB[1, t, s]) # Breeders to nonbreeders
        # Abundance of breeders
        N[1, 5, t+1, s] ~ dbin(mn.phiFY[1, t, s]*mn.psiFYB[1, t, s], NFY[1, t, s]) # Nestlings to second year breeders
        N[1, 6, t+1, s] ~ dbin(mn.phiA[1, t, s]*mn.psiAB[1, t, s], NF[1, t, s]) # Nonbreeder to breeder
        N[1, 7, t+1, s] ~ dbin(mn.phiB[1, t, s]*(1-mn.psiBA[1, t, s]), NB[1, t, s]) # Breeder to breeder
        # Assign estimates for years with data to 
        # all scenarios
        for (v in 1:7){ 
        for (sc in 2:SC){
          N[sc, v, t+1, s] <- N[1, v, t+1, s]
        } } # sc
        }} # s t

  # Future projections- population model 
    for (sc in 1:SC){
      for (t in nyr:(nyr+K-1)){
        for (s in 1:nsite){
          # Number of wild born juvs
          N[sc, 1, t+1, s] ~ dpois( (NFY[sc, t, s]*mn.phiFY[sc, t, s]*mn.psiFYB[sc, t, s] + # first year breeders
                                       NF[sc, t, s]*mn.phiA[sc, t, s]*mn.psiAB[sc, t, s] + # nonbreeders to breeders
                                       NB[sc, t, s]*mn.phiB[sc, t, s]*(1-mn.psiBA[sc, t, s])) # breeders remaining
                                    *(mn.prod[sc, t+1, s]/2) ) # end Poisson
          # Abundance of nonbreeders
          ## Second year nonbreeders
          N[sc, 2, t+1, s] ~ dbin(mn.phiFY[sc, t, s]*(1-mn.psiFYB[sc, t, s]), NFY[sc, t, s]) # Nestlings to second year nonbreeders
          ## Adult nonbreeders
          N[sc, 3, t+1, s] ~ dbin(mn.phiA[sc, t, s]*(1-mn.psiAB[sc, t, s]), NF[sc, t, s]) # Nonbreeders to nonbreeders
          N[sc, 4, t+1, s] ~ dbin(mn.phiB[sc, t, s]*mn.psiBA[sc, t, s], NB[sc, t, s]) # Breeders to nonbreeders
          # Abundance of breeders
          ## Second year breeders
          N[sc, 5, t+1, s] ~ dbin(mn.phiFY[sc, t, s]*mn.psiFYB[sc, t, s], NFY[sc, t, s]) # Nestlings to second year breeders
          ## Adult breeders
          N[sc, 6, t+1, s] ~ dbin(mn.phiA[sc, t, s]*mn.psiAB[sc, t, s], NF[sc, t, s]) # Nonbreeder to breeder
          N[sc, 7, t+1, s] ~ dbin(mn.phiB[sc, t, s]*(1-mn.psiBA[sc, t, s]), NB[sc, t, s]) # Breeder to breeder
        }} # s t
    } # sc
    
    # Count likelihoods, state-space model, and observation process 
    for (t in 1:nyr){
      for (s in 1:nsite){
        constraint_data[t, s] ~ dconstraint( (N[1, 1, t, s] + hacked.counts[1, t, s]) >= 0 ) # Transfers translocated first-year females
        NFY[1, t, s] <- N[1, 1, t, s] + hacked.counts[1, t, s] # Transfers translocated first-year females
        NF[1, t, s] <- sum(N[1, 2:4, t, s]) # number of adult nonbreeder females
        NB[1, t, s] <- sum(N[1, 5:7, t, s]) # number of adult breeder females
        NAD[1, t, s] <- NF[1, t, s] + NB[1, t, s] # number of adults
        Ntot[1, t, s] <- sum(N[1, 1:7, t, s]) # total number of females
        countsAdults[t, s] ~ dpois(lamAD[1, t, s]) # adult females
        # constrain N1+hacked.counts to be >=0
        log(lamAD[1, t, s]) <- log(NAD[1, t, s]) + deltas[1]*effort2[t, s] + deltas[2]*effort2[t, s]^2
        log(lamFY[1, t, s]) <- log(N[1, 1, t, s]) + deltas[3]*effort2[t, s] + deltas[4]*effort2[t, s]^2
        for (sc in 2:SC){
          lamAD[sc, t, s] <- lamAD[1, t, s]
          lamFY[sc, t, s] <- lamFY[1, t, s]
          NFY[sc, t, s] <- NFY[1, t, s]
          NF[sc, t, s] <- NF[1, t, s]
          NB[sc, t, s] <- NB[1, t, s]
          NAD[sc, t, s] <- NAD[1, t, s]
          Ntot[sc, t, s] <- Ntot[1, t, s]
        } # sc
      }# s
      # First-years at different sites have different distributions
      # for better model fit
      countsFY[t, 1] ~ dpois(lamFY[1, t, 1]) # doesn't have any zeroes so poisson
      countsFY[t, 2] ~ dnegbin(pp[t], r) # first year females, includes translocated/hacked
      pp[t] <- r/(r+(lamFY[1, t, 2] ))
    } # t
    r ~ dexp(0.05)
    
   # Future projections of counts 
    for(sc in 1:SC){
      for (t in (nyr+1):(nyr+K)){
        for (s in 1:nsite){
          # constrain N1+hacked.counts to be >=0, and allow it to shrink as the population shrinks
          NFY[sc, t, s] <- step( (N[sc, 1, t, s] + hacked.counts[sc, t, s]) ) * (N[sc, 1, t, s] + hacked.counts[sc, t, s]) 
          NF[sc, t, s] <- sum(N[sc, 2:4, t, s]) # number of adult nonbreeder females
          NB[sc, t, s] <- sum(N[sc, 5:7, t, s]) # number of adult breeder females
          NAD[sc, t, s] <- NF[sc, t, s] + NB[sc, t, s] # number of adults
          Ntot[sc, t, s] <- sum(N[sc, 1:7, t, s]) # total number of females
        }# s
      }  }# t sc
    
    ################################
    # Likelihood for survival
    ################################ 
    # Calculate yearly averages for integration
    for (t in 1:(nyr-1)){
      for (s in 1:nsite){
        for (xxxx in 1:surv.end[t,s]){
          # need to reorder because nimble doesn't 
          # handle nonconsecutive indices
          # yrind.surv is a matrix of indices for each site
          phiFY2[ t, s, xxxx] <- phiFY[ yrind.surv[xxxx,t,s], t]
          phiA2[ t, s, xxxx] <- phiA[ yrind.surv[xxxx,t,s], t]
          phiB2[ t, s, xxxx] <- phiB[ yrind.surv[xxxx,t,s], t]
          psiFYB2[ t, s, xxxx] <- psiFYB[ yrind.surv[xxxx,t,s], t]
          psiAB2[ t, s, xxxx] <- psiAB[ yrind.surv[xxxx,t,s], t]
          psiBA2[ t, s, xxxx] <- psiBA[ yrind.surv[xxxx,t,s], t]
          pA2[ t, s, xxxx] <- pA[ yrind.surv[xxxx,t,s], t]
          pB2[ t, s, xxxx] <- pB[ yrind.surv[xxxx,t,s], t]
        } # xxxx
        mn.phiFY[1, t, s] <- mean( phiFY2[ t, s, 1:surv.end[t,s] ] ) 
        mn.phiA[1, t, s] <- mean( phiA2[ t, s, 1:surv.end[t,s] ] )
        mn.phiB[1, t, s] <- mean( phiB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiFYB[1, t, s] <- mean( psiFYB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiAB[1, t, s] <- mean( psiAB2[ t, s, 1:surv.end[t,s] ] )
        mn.psiBA[1, t, s] <- mean( psiBA2[ t, s, 1:surv.end[t,s] ] )
        mn.pA[1, t, s] <- mean( pA2[ t, s, 1:surv.end[t,s] ] )
        mn.pB[1, t, s] <- mean( pB2[ t, s, 1:surv.end[t,s] ] )
        for (sc in 2:SC){
          mn.phiFY[sc, t, s] <- mn.phiFY[1, t, s] 
          mn.phiA[sc, t, s] <- mn.phiA[1, t, s] 
          mn.phiB[sc, t, s] <- mn.phiB[1, t, s] 
          mn.psiFYB[sc, t, s] <- mn.psiFYB[1, t, s]
          mn.psiAB[sc, t, s] <- mn.psiAB[1, t, s]
          mn.psiBA[sc, t, s] <- mn.psiBA[1, t, s]
          mn.pA[sc, t, s] <- mn.pA[1, t, s]
          mn.pB[sc, t, s] <- mn.pB[1, t, s]
        }
      }}
    
    for (i in 1:nind){
      for (t in 1:(nyr-1)){
        #Survival
        logit(phiFY[i,t]) <- eta[1, site[i,t],t] + eps[1,t] + 
          lmus[1, site[i,t]] + betas[1]*hacked[i]  # first year
        logit(phiA[i,t]) <- eta[2, site[i,t],t] + eps[2,t] + 
          lmus[2, site[i,t]] #+  betas[2]*hacked[i] # nonbreeder
        logit(phiB[i,t]) <- eta[3, site[i,t],t] + eps[3,t] + 
          lmus[3, site[i,t]] + betas[3]*hacked[i] # breeder
        #Recruitment
        logit(psiFYB[i,t]) <- eta[4, site[i,t],t] + eps[4,t] + 
          lmus[4, site[i,t]] #+ betas[4]*hacked[i] # first year to breeder
        logit(psiAB[i,t]) <- eta[5, site[i,t],t] + eps[5,t] + 
          lmus[5, site[i,t]] + betas[5]*hacked[i] # nonbreeder to breeder
        logit(psiBA[i,t]) <- #eta[6, site[i,t],t] + eps[6,t] + 
          lmus[6, site[i,t]] #+ betas[6]*hacked[i] # breeder to nonbreeder
        #Re-encounter
        logit(pA[i,t]) <- eta[7, site[i,t],t] + eps[7,t] + 
          lmus[7, site[i,t]] + betas[7]*hacked[i] +
          deltas[5]*effort2[t, site[i,t]] + deltas[6]*effort2[t, site[i,t]]^2 # resight of nonbreeders
        logit(pB[i,t]) <- eta[8, site[i,t],t] + eps[8,t] + 
          lmus[8, site[i,t]] + #betas[8]*hacked[i] + 
          deltas[7]*effort2[t, site[i,t]] + deltas[8]*effort2[t, site[i,t]]^2 # resight of breeders
      }#t
    }#i
    
    # Future projections of survival, recruitment, and detection
    for(sc in 1:SC){ 
      for (t in nyr:(nyr+K)){
        for (s in 1:nsite){
          # Calculate percent hacked
          # assigns a one (1.0 or 100%) if abundance is less than 10
          p.hackedFY[sc, t, s] <- step( NFY[sc, t, s]-10 ) * (10+0.001)/(NFY[sc, t, s]+0.001) + # NB>=10, calculate proportion
                                  step( NFY[sc, t, s]-1 ) * 1-step( NFY[sc, t, s]-10 )  # NB>1 and NB<10 100%
          p.hackedA[sc, t, s] <- step( NF[sc, t, s]-10 ) * (10+0.001)/(NF[sc, t, s]+0.001) + # NB>=10, calculate proportion
                                  step( NF[sc, t, s]-1 ) * 1-step( NF[sc, t, s]-10 )  # NB>1 and NB<10 100%
          p.hackedB[sc, t, s] <- step( NB[sc, t, s]-10 ) * (10+0.001)/(NB[sc, t, s]+0.001) + # NB>=10, calculate proportion
                                  step( NB[sc, t, s]-1 ) * 1-step( NB[sc, t, s]-10 )  # NB>1 and NB<10 100%
        # Hacking/translocation only affects birds
        # moved to Punta Cana (site 2) hence- equals(s, 2)
        #Survival
        logit(mn.phiFY[sc, t, s]) <- eta[1, s, t] + eps[1, t] + 
          lmus[1, s] + betas[1]*hacked2[sc]*p.hackedFY[sc, t, s]*equals(s, 2)  # first year
        logit(mn.phiA[sc, t, s]) <- eta[2, s, t] + eps[2, t] + 
          lmus[2, s] #+  betas[2]*hacked2[sc]*p.hackedA[sc, t, s]*equals(s, 2) # nonbreeder
        logit(mn.phiB[sc, t, s]) <- eta[3, s, t] + eps[3, t] + 
          lmus[3, s] + betas[3]*hacked2[sc]*p.hackedB[sc, t, s]*equals(s, 2) # breeder
        #Recruitment
        logit(mn.psiFYB[sc, t, s]) <- eta[4, s, t] + eps[4, t] + 
          lmus[4, s] #+ betas[4]*hacked2[sc]*p.hackedFY[sc, t, s]*equals(s, 2) # first year to breeder
        logit(mn.psiAB[sc, t, s]) <- eta[5, s, t] + eps[5, t] + 
          lmus[5, s] + betas[5]*hacked2[sc]*p.hackedA[sc, t, s]*equals(s, 2) # nonbreeder to breeder
        logit(mn.psiBA[sc, t, s]) <- #eta[6, site[i,t],t] + eps[6,t] + 
          lmus[6, s] #+ betas[6]*hacked[i] # breeder to nonbreeder
        #Re-encounter
        logit(mn.pA[sc, t, s]) <- eta[7, s, t] + eps[7, t] + 
          lmus[7, s] + betas[7]*hacked2[sc]*p.hackedA[sc, t, s]*equals(s, 2) # resight of nonbreeders
        logit(mn.pB[sc, t, s]) <- eta[8, s, t] + eps[8, t] + 
          lmus[8, s] #+ betas[8]*hacked2[sc]*p.hackedB[sc, t, s]*equals(s, 2) # resight of breeders
        } # s
      } # t
    } # sc
    
    # Define state-transition and observation matrices
    for (i in 1:nind){  
      # Define probabilities of state S(t+1) given S(t)
      for (t in first[i]:(nyr-1)){
        ps[1,i,t,1] <- 0
        ps[1,i,t,2] <- phiFY[i,t] * (1-psiFYB[i,t]) 
        ps[1,i,t,3] <- phiFY[i,t] * psiFYB[i,t]
        ps[1,i,t,4] <- (1-phiFY[i,t])
        
        ps[2,i,t,1] <- 0
        ps[2,i,t,2] <- phiA[i,t] * (1-psiAB[i,t])
        ps[2,i,t,3] <- phiA[i,t] * psiAB[i,t]
        ps[2,i,t,4] <- (1-phiA[i,t])
        
        ps[3,i,t,1] <- 0
        ps[3,i,t,2] <- phiB[i,t] * psiBA[i,t]
        ps[3,i,t,3] <- phiB[i,t] * (1-psiBA[i,t])
        ps[3,i,t,4] <- (1-phiB[i,t])
        
        ps[4,i,t,1] <- 0
        ps[4,i,t,2] <- 0
        ps[4,i,t,3] <- 0
        ps[4,i,t,4] <- 1
        
        # Define probabilities of O(t) given S(t)
        po[1,i,t,1] <- 1 
        po[1,i,t,2] <- 0
        po[1,i,t,3] <- 0
        po[1,i,t,4] <- 0
        
        po[2,i,t,1] <- 0
        po[2,i,t,2] <- pA[i,t]
        po[2,i,t,3] <- 0
        po[2,i,t,4] <- (1-pA[i,t])
        
        po[3,i,t,1] <- 0
        po[3,i,t,2] <- 0
        po[3,i,t,3] <- pB[i,t]
        po[3,i,t,4] <- (1-pB[i,t])
        
        po[4,i,t,1] <- 0
        po[4,i,t,2] <- 0
        po[4,i,t,3] <- 0
        po[4,i,t,4] <- 1
      } #t
    } #i
    
    # Likelihood 
    for (i in 1:nind){
      # Define latent state at first capture
      z[i,first[i]] <- y[i,first[i]]
      for (t in (first[i]+1):nyr){
        # State process: draw S(t) given S(t-1)
        z[i,t] ~ dcat(ps[z[i,t-1], i, t-1, 1:4])
        # Observation process: draw O(t) given S(t)
        y[i,t] ~ dcat(po[z[i,t], i, t-1, 1:4])
      } #t
    } #i
  } )

#**********************
#* Function to run model in NIMBLE
#**********************
run_pva <- function(info, datl, constl, code){
  library('nimble')
  library('coda')
  # function for multivariate normal
  uppertri_mult_diag <- nimbleFunction(
    run = function(mat = double(2), vec = double(1)) {
      returnType(double(2))
      p <- length(vec)
      out <- matrix(nrow = p, ncol = p, init = FALSE)
      for(i in 1:p){
        out[ ,i] <- mat[ ,i] * vec[i]
      }
      return(out)
    })
  assign('uppertri_mult_diag', uppertri_mult_diag, envir = .GlobalEnv)
  
  params <- c(
    # pop growth 
    "lambda",
    # prod
    "lmu.prod", "gamma", "rr", "mn.prod", 
    # survival 
    "mus", "lmus", "betas", "deltas",
    # abundance
    "NB", "NF", "NFY", "N", "NAD", "Ntot",
    "r",
    # error terms
    "eps", "sds", "Ustar", "U", "R",
    "eta", "sds2", "Ustar2", "U2", "R2",
    # yearly summaries
    'mn.phiFY', 'mn.phiA', 'mn.phiB',
    'mn.psiFYB', 'mn.psiAB', 'mn.psiBA',
    'mn.pA', 'mn.pB',
    # pva
    "extinct", "extinctAD", "extinctB"
  )
  #n.chains=1; n.thin=1; n.iter=500; n.burnin=100
  n.chains=1; n.thin=200; n.iter=600000; n.burnin=400000
  
  mod <- nimbleModel(code, 
                     constants = constl, 
                     data = datl, 
                     inits = info$inits, 
                     buildDerivs = FALSE,
                     calculate=T 
  ) 
  
  cmod <- compileNimble(mod, showCompilerOutput = TRUE)
  confhmc <- configureMCMC(mod)
  confhmc$setMonitors(params)
  hmc <- buildMCMC(confhmc)
  chmc <- compileNimble(hmc, project = mod, 
                        resetFunctions = TRUE,
                        showCompilerOutput = TRUE)
  
  post <- runMCMC(chmc,
                  niter = n.iter, 
                  nburnin = n.burnin,
                  nchains = n.chains,
                  thin = n.thin,
                  samplesAsCodaMCMC = T, 
                  setSeed = info$seed)
  
  return(post)
} # run_pva function end

#*****************
#* Run chains in parallel
#*****************
this_cluster <- makeCluster(cpus)
post <- parLapply(cl = this_cluster, 
                  X = par_info_pva[1:cpus], 
                  fun = run_pva, 
                  datl = datl, 
                  constl = constl, 
                  code = mycode)
stopCluster(this_cluster)
save(post, mycode, seeds, cpus,
     file="/bsuscratch/brianrolek/riha_ipm/outputs/pva.rdata")

2. Figures and tables

2.1 Plot model estimates

library ('MCMCvis')
library ('coda')
library ('ggplot2')
library('reshape2')
library ("tidybayes")
library ('bayestestR')
library ("ggpubr")
library("viridis")
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm_longrun.rdata")
load("data/data.rdata")
out <- lapply(post, as.mcmc) 
outp <- MCMCpstr(out, type="chains")

#********************
#* Calculate summary stats
#********************
# Population growth rates
mn.lambda <- apply(outp$lambda, c(2,3), mean)
(md.lambda <- apply(mn.lambda, 1, median))
## [1] 1.000060 1.304276
(hdi95.lambda <- apply(mn.lambda, 1, HDInterval::hdi))
##            [,1]     [,2]
## lower 0.9479604 1.209558
## upper 1.0348278 1.427075
# Overall averages
mn.mus <- apply(outp$mus, c(1,3), mean)
(md.mus <- apply(mn.mus, 1, median))
## [1] 0.33112222 0.93945299 0.88785788 0.11015826 0.40630345 0.07551559 0.20468205
## [8] 0.78345859
(hdi.mus <- apply(mn.mus, 1, HDInterval::hdi))
##            [,1]      [,2]      [,3]       [,4]      [,5]      [,6]       [,7]
## lower 0.1721033 0.8122186 0.7967305 0.01096642 0.2414542 0.0327038 0.04307028
## upper 0.5341632 0.9994355 0.9647374 0.30829937 0.5460728 0.1214192 0.42815018
##            [,8]
## lower 0.5467521
## upper 0.9503702
# Compare between sites
df.sites <- data.frame(mus=1:8, pd=rep(NA,8), greater=NA)
for (i in 1:8){
  first <- median(outp$mus[i, 1, ]) > median(outp$mus[i, 2, ])
  df.sites$greater[i] <- ifelse(first, "LH", "PC")
  if (first){
    diff <- outp$mus[i, 1, ]-outp$mus[i, 2, ]
    df.sites$pd[i] <- mean(diff>0) |> round(2) 
  } else{
    diff <- outp$mus[i, 2, ]-outp$mus[i, 1, ]
    df.sites$pd[i] <- mean(diff>0) |> round(2)
  }
}
df.sites
(md.mus <- apply(outp$mus, c(1,2), median)) |> round(2)
##      [,1] [,2]
## [1,] 0.39 0.26
## [2,] 0.97 0.93
## [3,] 0.91 0.87
## [4,] 0.12 0.07
## [5,] 0.11 0.70
## [6,] 0.12 0.03
## [7,] 0.06 0.34
## [8,] 0.71 0.89
(hdi.mus <- apply(outp$mus, c(1,2), HDInterval::hdi)) |> round(2)
## , , 1
## 
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## lower 0.19 0.87 0.83 0.01 0.04 0.05  0.0 0.34
## upper 0.65 1.00 0.98 0.35 0.19 0.20  0.2 0.96
## 
## , , 2
## 
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## lower 0.07 0.68 0.74  0.0 0.41 0.00 0.04  0.6
## upper 0.52 1.00 0.97  0.4 0.95 0.08 0.74  1.0
#*******************
#* Plot IPM results
#*******************
#Plot population size
Ntot <- outp$Ntot[1:13,,] # thin by 10 to speed up plotting
lN <- melt(Ntot)
colnames(lN) <- c("Time", "Site.ind", "Iter", "Abundance")
lN$Year <- lN$Time + 2010
lN$Site <- ifelse(lN$Site.ind==1, "Los Haitises", "Punta Cana")

med_hdis <- function(x, cm){
            df <- data.frame(y=median(x),
                             ymin=HDInterval::hdi(x, credMass=cm)[1],
                             ymax=HDInterval::hdi(x, credMass=cm)[2] )
}

counts.tot <- datl$countsAdults+datl$countsFY+constl$hacked.counts[,-3]
df.counts <- melt(counts.tot)
colnames(df.counts) <- c("Year", "Site_ab", "Count")
df.counts$Site <- ifelse(df.counts$Site_ab=="LHNP", "Los Haitises", "Punta Cana")

p1 <- ggplot() + theme_minimal(base_size=14) + 
  geom_line(data=lN, aes(x=Year, y=Abundance, group=Iter), 
            color="gray30", linewidth=0.1, alpha=0.01 ) +
  stat_summary(data=lN, aes(x=Year, y=Abundance), 
               geom="errorbar", width=0, 
               fun.data =med_hdis, fun.args=list(cm=0.95) ) +
  stat_summary(data=lN, aes(x=Year, y=Abundance), 
               geom="errorbar", width=0, linewidth=1,
               fun.data =med_hdis, fun.args=list(cm=0.80) ) +
  stat_summary(data=lN, aes(x=Year, y=Abundance), 
               geom="line", linewidth=1,
               fun.data =function(x){data.frame(y=median(x))} ) +
  geom_line(data=df.counts, aes(x=Year, y=Count), 
            col="black", linewidth=1, linetype=2 ) +
  facet_wrap("Site", scales="free_y") +
  xlab("Year") + ylab("Number")
p1

# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Abundance and counts.tiff",
#        plot=p1,
#        device="tiff",
#        width=6, height=4, dpi=300)
# Plot life stage structure
mdFY <-  apply(outp$NFY[1:13,,], c(1,2), median) 
mdB <-  apply(outp$NB[1:13,,], c(1,2), median) 
mdF <-  apply(outp$NF[1:13,,], c(1,2), median) 
lFY <- melt(mdFY)
lB <- melt(mdB)
lF <- melt(mdF)
lFY$Stage <- "First-year"
lB$Stage <- "Breeder"
lF$Stage <- "Nonbreeder"
ldat <- rbind(lFY, lB, lF)
colnames(ldat)[1:3] <- c("Year", "Sitenum", "Number") 
ldat$Site <- ifelse(ldat$Sitenum==1, "Los Haitises", "Punta Cana")
ldat$year <- ldat$Year+2010

# Use median number of females in each stage
# to plot an approximate population structure
p3 <- ggplot(ldat, aes(fill=Stage, y=as.numeric(Number), x=year)) + 
  geom_bar(position="fill", stat="identity") +
  ylab("Proportion of population") + xlab("Year") +
  facet_wrap("Site") + 
  theme_minimal(base_size=14)+ theme(legend.position="top") +
  scale_fill_viridis_d(option = "mako")

p4 <- ggplot(ldat, aes(fill=Stage, y=as.numeric(Number), x=year)) + 
  geom_bar(position="stack", stat="identity") +
  ylab("Median number of females") + xlab("Year") +
  facet_wrap("Site", scales = "free_y") + 
  theme_minimal(base_size=14 ) + theme(legend.position="top") + 
  scale_fill_viridis_d(option = "mako") 

# plot with uncertainty
lFY <- melt(outp$NFY)
lB <- melt(outp$NB)
lF <- melt(outp$NF)
colnames(lFY) <- colnames(lB) <- colnames(lF) <- c("Time", "Site", "Iteration", "Abundance")
lFY$Year <- lFY$Time + 2010
lB$Year <- lB$Time + 2010
lF$Year <- lF$Time + 2010
lFY$Stage <- "First-year"
lB$Stage <- "Breeder"
lF$Stage <- "Nonbreeder"
lALL <- rbind(lFY, lB, lF)
lALL$Site <- ifelse(lALL$Site==1, "Los Haitises", "Punta Cana")

mdFY <-  apply(outp$NFY[1:13,,], c(1,2), median) 
mdB <-  apply(outp$NB[1:13,,], c(1,2), median) 
mdF <-  apply(outp$NF[1:13,,], c(1,2), median) 
hdiFY <-  apply(outp$NFY[1:13,,], c(1,2), HDInterval::hdi) 
hdiB <-  apply(outp$NB[1:13,,], c(1,2), HDInterval::hdi) 
hdiF <-  apply(outp$NF[1:13,,], c(1,2), HDInterval::hdi) 
medFY <- melt(mdFY)
medB <- melt(mdB)
medF <- melt(mdF)
hdisFY <- melt(hdiFY)
hdisB <- melt(hdiB)
hdisF <- melt(hdiF)
dfstages <- rbind(medFY, medB, medF)
dfhdis <- rbind(hdisFY, hdisB, hdisF)
dfstages$LHDI <- dfhdis[dfhdis$Var1=="lower",]
dfstages$UHDI <- dfhdis[dfhdis$Var1=="upper",]
cols <- mako(3)

med_hdis <- function(x, cm){
  df <- data.frame(y=median(x),
                   ymin=HDInterval::hdi(x, credMass=cm)[1],
                   ymax=HDInterval::hdi(x, credMass=cm)[2] )
}

p5 <- ggplot() +  theme_minimal(base_size=14) +
  theme(legend.position="none") +
  geom_line(data=lALL, aes(x=Year, y=Abundance, 
                           group=Iteration, color=Stage,
                           alpha=Stage), 
            linewidth=0.1) +
  stat_summary(data=lALL, aes(x=Year, y=Abundance),
               geom="errorbar", width=0,
               fun.data =med_hdis, fun.args=list(cm=0.95) ) +
  stat_summary(data=lALL, aes(x=Year, y=Abundance),
               geom="errorbar", width=0, linewidth=1,
               fun.data =med_hdis, fun.args=list(cm=0.80) ) +
  stat_summary(data=lALL, aes(x=Year, y=Abundance),
               geom="line", linewidth=1,
               fun.data =function(x){data.frame(y=median(x))} ) +
  scale_color_manual( values = c("First-year"=cols[2], "Breeder"=cols[1], "Nonbreeder"=cols[3])  ) +
  scale_alpha_manual( values = c("First-year"=0.015, "Breeder"=0.0075, "Nonbreeder"=0.15) ) +
  facet_wrap(Stage~Site, scales="free_y", 
             shrink=TRUE, ncol=2) +
  xlab("Year") + ylab("Number of females") +
  theme(
    strip.background = element_blank(),
    strip.text.x = element_blank()
  )

p35 <- ggarrange(p3, p5, ncol=1, nrow=2)
p35

# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Stage structure_abundance.tiff",
#        plot=p35,
#        device="tiff",
#        width=6, height=8, dpi=300)
# plot survival, recruitment, and detection
mus.mat <- outp$mus
dimnames(mus.mat)[[1]] <- c("FY", 
                            "A\nSurvival", "B", 
                            "FY:B",
                            "A:B\nRecruitment", 
                            "B:A",
                            "A", "B\nDetection")
dimnames(mus.mat)[[2]] <- c("Los Haitises", "Punta Cana")
lmus <- melt(mus.mat)
colnames(lmus) <- c("Parameter", "Site", "Iter", "value")

p6 <- ggplot(data= lmus, aes(x=value, y=Parameter )) +
  stat_halfeye(point_interval=median_hdi, .width = c(.95, .80)) +
  coord_flip() +
  facet_wrap("Site") +
  xlab("Probability") + ylab("Parameter") +
  theme_bw(base_size=14) +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank()) 

p6

# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Survival Recruitment Detection.tiff",
#        plot=p5,
#        device="tiff",
#        width=8, height=4, dpi=300)


# Plot predicted survival, recruitment, and detection
# with effects from translocation and hacking
# Birds were only hacked from LHNP to PC here
# so we only predict values for PC
# pred.mus <- array(NA, dim=dim(outp$mus))
# for (m in 1:8){
#   for (tr in 1:2){
#     pred.mus[m,tr,] <- plogis( outp$lmus[m,2,] + outp$betas[m,]*c(0,1)[tr] ) 
#   }}
# 
# dimnames(pred.mus) <- list(param=c("FY", "A", "B", "FY:B", "A:B", "B:A", "A", "B"),
#                            trans=c("Not translocated", "Translocated"), 
#                            iter=1:10000)
# lm <- melt(pred.mus)
# 
# p5 <- ggplot(data= lm, aes(x=value, y=param )) +
#   stat_halfeye(point_interval=median_hdi, .width = c(.95, .80)) +
#   coord_flip() +
#   facet_wrap("trans") +
#   xlab("Probability") + ylab("Parameter") +
#   theme_bw(base_size=14) +
#   theme(panel.grid.major.x = element_blank(),
#         panel.grid.minor.x = element_blank()) 
betas.pd <- apply(outp$betas, 1, pd)
ind <- which (betas.pd>0.975) # which are significant
pred.mus <- array(NA, dim=c(dim(outp$mus)))
for (m in ind){
  for (tr in 1:2){
    pred.mus[m,tr,] <- plogis( outp$lmus[m,2,] + outp$betas[m,]*c(0,1)[tr] ) 
  }}
# treatment, mu, site, iter
mus.md <- apply(pred.mus, c(1,2), median)
mus.HDI80 <- apply(pred.mus, c(1,2), HDInterval::hdi, credMass=0.8)
mus.HDI95 <- apply(pred.mus, c(1,2), HDInterval::hdi)
# tiff(height=4, width=6, units="in", res=300,
#      filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Translocation effects.tiff")
par(mar=c(3,4,4,1), mfrow=c(1,1))
for (tr in 1:2){
  if (tr==1){
    plot(1:4, c(0,0,1,1), 
         type="n",
         pch=c(16, 17)[tr], cex=2,
         ylim=c(0,1),
         xlim=c(0.5, 4.5), 
         cex.axis=1, cex.lab=1.25,
         ylab="Probability", xlab="",
         main="",
         xaxt="n", yaxt="n")
  abline(h= seq(0,1,by=0.1), col="gray90")
  abline(v= seq(0.5,6,by=1), col="gray90")
  points(1:4+c(-0.1,0.1)[tr], mus.md[ind,tr], 
         pch=c(16, 17)[tr], cex=1.5)
  }else{
    points(1:4+c(-0.1,0.1)[tr], mus.md[ind,tr], 
           pch=c(16, 15)[tr], cex=1.5)
  }}

axis(1, at=1:4, cex.axis=0.85,
     labels=c("First-year\nSurvival", 
              "Breeder\nSurvival", 
              "Nonbreeder\n to Breeder\nRecruitment", 
              "Nonbreeder\nDetection"), las=1, padj=0.5)
axis(2, at=seq(0,1, by=0.1), cex.axis=1,
     labels=c(0, NA,NA,NA,NA,0.5, NA,NA,NA,NA, 1))
for (m in 1:length(ind)){
  for (tr in 1:2){
    lines(rep(c(1:4)[m] + c(-0.1,0.1)[tr],2), mus.HDI95[1:2,ind[m],tr], lwd=2)
    lines(rep(c(1:4)[m] + c(-0.1,0.1)[tr],2), mus.HDI80[1:2,ind[m],tr], lwd=4)
  }}
legend(x=2.5,y=1.25, pch=c(16,15), pt.cex=1.5, cex=1, xpd=NA, horiz=T, 
       legend=c("Not translocated", "Translocated" ),
       xjust=0.5)

#dev.off()
# Plot fecundity in response to nest treatments
f <- exp(outp$lmu.prod)
f.md <- apply(f, 1, median)
f.HDI80 <- apply(f, 1, HDInterval::hdi, credMass=0.8)
f.HDI95 <- apply(f, 1, HDInterval::hdi)

# Calculate treatment effects on fecundity
f.pred <- array(NA, dim=dim(outp$lmu.prod))
for (s in 1:2){
  f.pred[s,] <- exp(outp$lmu.prod[s,] + outp$gamma[,1])
} # s
f2.md <- apply(f.pred, 1, median)
f2.HDI80 <- apply(f.pred, 1, HDInterval::hdi, credMass=0.8)
f2.HDI95 <- apply(f.pred, 1, HDInterval::hdi)

# tiff(height=4, width=6, units="in", res=300,
#      filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Fecundity_treatment.tiff")
par(mar=c(3,5,3,1), mfrow=c(1,1))
plot(c(1, 2)-0.1, f.md, 
     ylim=c(min(f.HDI95), max(f2.HDI95)),
     xlim=c(0.5, 2.5), 
     ylab="Productivity", xlab="",
     main="",
     cex.axis=1.5, cex.lab=1.5,
     xaxt="n", yaxt="n", type="n")
points(c(1, 2)-0.1, f.md, pch=16, cex=1.5)
abline(h=seq(0, 1.6, by=0.1), col="gray90")
abline(v=1.5, col="black")
axis(1, at=c(1,2), cex.axis=1,
     labels=c("Los Haitises\nNational Park","Punta Cana"),
     padj=0.5)
axis(2, at=seq(0, 1.4, by=0.1), 
     cex.axis=1,
     labels=c(0, NA, NA, NA, 0.4, NA, NA, NA, 0.8, NA, NA, NA, 1.2, NA, NA))
lines(c(1,1)-0.1, f.HDI95[,1], lwd=2)
lines(c(2,2)-0.1, f.HDI95[,2], lwd=2)
lines(c(1,1)-0.1, f.HDI80[,1], lwd=4)
lines(c(2,2)-0.1, f.HDI80[,2], lwd=4)

points(c(1.1, 2.1), f2.md, 
       pch=18, cex=2.5)

lines(c(1,1) +0.1, f2.HDI95[,1], lwd=3)
lines(c(2,2) +0.1, f2.HDI95[,2], lwd=3)
lines(c(1,1) +0.1, f2.HDI80[,1], lwd=6)
lines(c(2,2) +0.1, f2.HDI80[,2], lwd=6)
legend(x=1.5,y=1.65, pch=c(16,18), pt.cex=c(1.5,2.5), cex=1,
       legend=c("Untreated", "Treated" ), 
       horiz=TRUE, xjust=0.5, xpd=NA)

#dev.off()

# Calculate probability of direction
f.diff <- f[1,]-f[2,]
pd(f.diff)
## Probability of Direction: 1.00
f.md 
## lmu.prod[1] lmu.prod[2] 
##   0.3659702   0.2201483
f.HDI95 
##       lmu.prod[1] lmu.prod[2]
## lower   0.2877536   0.1581736
## upper   0.4512852   0.2860523
f.diff2 <- f.pred[1,]-f.pred[2,]
pd(f.diff2)
## Probability of Direction: 1.00
f2.md 
## [1] 1.1107340 0.6681587
f2.HDI95
##            [,1]      [,2]
## lower 0.8733436 0.4800632
## upper 1.3696686 0.8681799
median(f.pred[1,]/f[1,])
## [1] 3.03504
median(f.pred[2,]/f[2,])
## [1] 3.03504
#*******************
# Correlations between demographic rates
#*******************
apply(outp$R, c(1,2), pd) 
##          [,1]    [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]    [,9]
##  [1,] 1.00000 0.51395 0.56615 0.61530 0.55105 0.50915 0.54825 0.50140 0.55795
##  [2,] 0.51395 1.00000 0.50495 0.54665 0.54755 0.51475 0.58975 0.50035 0.51885
##  [3,] 0.56615 0.50495 1.00000 0.58120 0.61845 0.51340 0.63790 0.51505 0.53490
##  [4,] 0.61530 0.54665 0.58120 1.00000 0.60075 0.50615 0.71585 0.56350 0.68935
##  [5,] 0.55105 0.54755 0.61845 0.60075 1.00000 0.50560 0.57150 0.54245 0.55390
##  [6,] 0.50915 0.51475 0.51340 0.50615 0.50560 1.00000 0.51845 0.50560 0.50450
##  [7,] 0.54825 0.58975 0.63790 0.71585 0.57150 0.51845 1.00000 0.52575 0.74405
##  [8,] 0.50140 0.50035 0.51505 0.56350 0.54245 0.50560 0.52575 1.00000 0.53540
##  [9,] 0.55795 0.51885 0.53490 0.68935 0.55390 0.50450 0.74405 0.53540 1.00000
apply(outp$R2, c(1,2), pd) 
##          [,1]    [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]    [,9]
##  [1,] 1.00000 0.51055 0.51580 0.52700 0.53725 0.50115 0.50535 0.52995 0.52010
##  [2,] 0.51055 1.00000 0.53865 0.50220 0.50155 0.50645 0.50730 0.56050 0.50395
##  [3,] 0.51580 0.53865 1.00000 0.53130 0.54490 0.51285 0.55790 0.56290 0.52050
##  [4,] 0.52700 0.50220 0.53130 1.00000 0.53680 0.50315 0.55635 0.65250 0.55840
##  [5,] 0.53725 0.50155 0.54490 0.53680 1.00000 0.50265 0.52515 0.52370 0.55790
##  [6,] 0.50115 0.50645 0.51285 0.50315 0.50265 1.00000 0.50840 0.51085 0.51740
##  [7,] 0.50535 0.50730 0.55790 0.55635 0.52515 0.50840 1.00000 0.73010 0.52495
##  [8,] 0.52995 0.56050 0.56290 0.65250 0.52370 0.51085 0.73010 1.00000 0.57520
##  [9,] 0.52010 0.50395 0.52050 0.55840 0.55790 0.51740 0.52495 0.57520 1.00000
median(outp$R2[7,8,]) 
## [1] -0.1869805
#*******************
# Correlations between demographics and population growth rates
#*******************
lam.m <- apply(outp$lambda[1:12,,], c(1,2), median)
lam.hdi <- apply(outp$lambda[1:12,,], c(1,2), HDInterval::hdi)
par(mfrow=c(1,2))
plot(2012:2023, lam.m[,1], type="b", pch=1, 
     ylab="Population growth rate", xlab="Year", 
     ylim=c(min(lam.hdi[,,1]), max(lam.hdi[,,1])),
     main="Los Haitises")
abline(h=1, lty=2)
segments(x0=2012:2023, x1=2012:2023, 
         y0 = lam.hdi[1,,1], y1= lam.hdi[2,,1])

plot(2012:2023, lam.m[,2], type="b", pch=1, lty=1,
     ylab="Population growth rate", xlab="Year",
     ylim=c(min(lam.hdi[,,2]), max(lam.hdi[,,2])),
     main="Punta Cana")
abline(h=1, lty=2)
segments(x0=2012:2023, x1=2012:2023, 
         y0 = lam.hdi[1,,2], y1= lam.hdi[2,,2])

# create a function to plot correlations
# between demographics and population growth rates
plot.cor <- function (lam.post, x.post, x.lab, ind.x=1:12, site=1, yaxt="b"){
  # calculate the correlation coefficient 
  # over each iteration to propagate error
  lambda <- lam.post[1:12, site,]
  x <- x.post[ind.x, site, ]
  df <- data.frame(ats1=c(0.8, 0.9, 1.0, 1.1, 1.2, 1.3),
                   ats2=c(1.0, 1.5, 2.0, 2.5, 3.0, 3.5),
                   labs1=c("0.8", NA, "1.0", NA, "1.2", NA),
                   labs2=c("1.0", NA, "2.0", NA, "3.0", NA)
  )
  if(site==1){ df <- df[, c(1,3)]} else{
    df <- df[, c(2,4)]
  }
  lwd <- 1
  
  cor.post <- array(NA, dim=dim(lambda)[-1])
    for (i in 1:dim(lambda)[[2]]){
      cor.post[i] <- cor(lambda[,i], x[,i])
    }
  cor.df <- data.frame(median= median(cor.post) |> round(digits=2),
                       ldi=HDInterval::hdi(cor.post)[1] |> round(digits=2),
                       hdi=HDInterval::hdi(cor.post)[2] |> round(digits=2),
                       pd = pd(cor.post) |> round(digits=2))
  
  lam.m <- apply(lambda, 1, median)
  lam.hdi <- apply(lambda, 1, HDInterval::hdi)
  x.m <- apply(x, 1, median)
  x.hdi <- apply(x, 1, HDInterval::hdi)
    x.lims <- c(min(x.hdi[,]), max(x.hdi[,]))
    y.lims <- c(min(lam.hdi[,]), max(lam.hdi[,]))
    if(yaxt=="n"){
    plot(x.m, lam.m[1:12], 
         xlim= x.lims,
         ylim= y.lims,
         type="n", 
         ylab="", xlab=x.lab, cex.lab=1.5, 
         main=c("", "")[site],
         yaxt=yaxt, xaxt="n")
      axis(2, at=df[,1], labels=c(rep(NA, 5)))
    } else {
      plot(x.m, lam.m[1:12], 
           xlim= x.lims,
           ylim= y.lims,
           type="n", 
           ylab="", xlab=x.lab, cex.lab=1.5,  
           main=c("", "")[site], 
           yaxt="n", xaxt="n")
      axis(2, at=df[,1], labels=df[,2], las=1, cex.axis=1.5)
    }
    axis(1, cex.axis=1.5)
    points(x.m, lam.m, pch=1)
    segments(x0=x.hdi[1,], x1=x.hdi[2,], 
             y0 = lam.m, y1= lam.m, lwd=lwd)
    segments(x0=x.m, x1=x.m, 
             y0 = lam.hdi[1,], y1= lam.hdi[2,], lwd=lwd)
    xs <- c(x.lims[1], x.lims[1]+(x.lims[2]-x.lims[1]*1.5))
    ys <- c( (y.lims[2]-y.lims[1])+y.lims[1], (y.lims[2]-y.lims[1])+y.lims[1], 
             (y.lims[2]-y.lims[1])*0.6+y.lims[1], (y.lims[2]-y.lims[1])*0.6+y.lims[1])
    polygon(x=c(xs, rev(xs)), y=ys, col=alpha("white", 0.7), border=NA )
    text(x = x.lims[1], y = (y.lims[2]-y.lims[1])*0.9+y.lims[1], 
         paste("r = ", cor.df$median, " (",cor.df$ldi,", ", cor.df$hdi, ")", sep=""), 
         pos = 4, font = 3, cex = 1.5)
    text(x = x.lims[1], y = (y.lims[2]-y.lims[1])*0.7+y.lims[1], paste("P(r>0) = ", cor.df$pd, sep=""), 
         pos = 4, font = 3, cex = 1.5)
    
  }

# Plor correlations between population grrowth rates
# and demographics. 
# "r" is a correlation coefficient and represents
# the magnitude of the correlation
# P(r>0) is the probability of direction (similar to p-values)
# that is, the probability that an effect exists
# tiff(height=8, width=8, units="in", res=300,
#      filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//PopGrowthand Demographics1.tiff")
# par(mfrow=c(4,3), mar=c(4,1,1,1), oma=c(0,5,3,0))
par(mfrow=c(2,3), mar=c(4,1,1,1), oma=c(0,5,3,0))
plot.cor(outp$lambda, outp$mn.prod, x.lab="Fecundity", ind.x=2:13)
plot.cor(outp$lambda, outp$mn.phiFY, x.lab="First-year Survival", yaxt="n")
plot.cor(outp$lambda, outp$mn.phiA, x.lab="Nonbreeder Survival", yaxt="n")
plot.cor(outp$lambda, outp$mn.phiB, x.lab="Breeder Survival")
plot.cor(outp$lambda, outp$mn.psiFYB, x.lab="First-year to Breeder", yaxt="n")
plot.cor(outp$lambda, outp$mn.psiAB, x.lab="Nonbreeder to Breeder", yaxt="n")
mtext("Population growth rate", side=2, outer=TRUE, line=2.5, cex=2)
mtext("Los Haitises", side=3, outer=TRUE, cex=2)

#dev.off()

# tiff(height=4, width=8, units="in", res=300,
#      filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//PopGrowthand Demographics2.tiff")
par(mfrow=c(2,3), mar=c(4,1,1,1), oma=c(0,5,3,0))
plot.cor(outp$lambda, outp$mn.prod, x.lab="Fecundity", ind.x=2:13, site=2)
plot.cor(outp$lambda, outp$mn.phiFY, x.lab="First-year Survival", yaxt="n", site=2)
plot.cor(outp$lambda, outp$mn.phiA, x.lab="Nonbreeder Survival", yaxt="n", site=2)
plot.cor(outp$lambda, outp$mn.phiB, x.lab="Breeder Survival", site=2)
plot.cor(outp$lambda, outp$mn.psiFYB, x.lab="First-year to Breeder", yaxt="n", site=2)
plot.cor(outp$lambda, outp$mn.psiAB, x.lab="Nonbreeder to Breeder", yaxt="n", site=2)
mtext("Population growth rate", side=2, outer=TRUE, line=2.5, cex=2)
mtext("Punta Cana", side=3, outer=TRUE, cex=2)

#dev.off()
#*******************
#* Plot PVA results
#*******************
# Abundance of females at Los Haitises
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\pva.rdata")
out <- lapply(post, as.mcmc) 
outp <- MCMCpstr(out, type="chains")
library (viridis)
cols <- viridis(3)
pe <- apply(outp$extinct, c(1,2,3), mean)

# tiff(height=3, width=8, units="in", res=300,
#      filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Quasiextinction.tiff")
par(mar=c(2, 2, 1, 1), oma=c(3,3,0,0))
layout(matrix(c(1,1,2,2,3), 1, 5, byrow = TRUE))

plot(2024:2123, pe[1, ,1], 
     type="n", col=cols[1], lwd=2,
     ylim=c(0,1), main="Los Haitises",
     ylab="Quasi-extinction probability", xlab="Year",
     xaxt="n")
axis(1, at=seq(2020, 2120, by=10), 
     labels=c(2020, NA, 2040, NA, 2060, NA, 2080, NA, 2100, NA, 2120))
abline(h=seq(0, 1, by=0.1), col="gray90")
abline(v=seq(2030, 2130, by=10), col="gray90")
lines(2024:2123, pe[1, ,1], col=cols[1], lwd=2)
lines(2024:2123, pe[2, ,1], col=cols[2], lwd=2)
lines(2024:2123, pe[3, ,1], col=cols[3], lwd=2)
lines(2024:2123, pe[4, ,1], col=cols[1], lwd=2, lty=2)
lines(2024:2123, pe[5, ,1], col=cols[2], lwd=2, lty=2)
lines(2024:2123, pe[6, ,1], col=cols[3], lwd=2, lty=2)

plot(2024:2123, pe[1, ,2], 
     type="n", col=cols[1], lwd=2,
     ylim=c(0,1), main="Punta Cana",
     ylab="", xlab="Year",
     xaxt="n")
axis(1, at=seq(2020, 2120, by=10), 
     labels=c(2020, NA, 2040, NA, 2060, NA, 2080, NA, 2100, NA, 2120))
abline(h=seq(0, 1, by=0.1), col="gray90")
abline(v=seq(2030, 2130, by=10), col="gray90")
lines(2024:2123, pe[1, ,2], col=cols[1], lwd=2)
lines(2024:2123, pe[2, ,2], col=cols[2], lwd=2)
lines(2024:2123, pe[3, ,2], col=cols[3], lwd=2)
lines(2024:2123, pe[4, ,2], col=cols[1], lwd=2, lty=2)
lines(2024:2123, pe[5, ,2], col=cols[2], lwd=2, lty=3)
lines(2024:2123, pe[6, ,2], col=cols[3], lwd=2, lty=4)
plot.new()
legend(x=-0.3, y=0.5, title="Scenarios",
       legend=c("1 trans=0, nests=0", "2 trans=0, nests=10", "3 trans=0, nests=all",
                "4 trans=10, nests=0", "5 trans=10, nests=10", "6 trans=10, nests=all"),
       xpd=NA, col=c(cols, cols), lty=c(1,1,1,2,2,2), lwd=2,
       xjust=0, yjust=0.5)
mtext("Future year", 1,  outer=TRUE, line=1, adj=0.39)
mtext("Quasi-extinction probability", 2, outer=TRUE, line=1)

# dev.off()

# Time to extinction
# conditional on extinction threshold (D=X)
# D <- 3
# qextinct <- outp$Ntot <=3
# ext <- outp$extinct[,100,,]==1
# t.extinction <- apply(qextinct[])

2.2 Plots of finer population segments and model diagnostics

load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm_longrun.rdata")
load("data/data.rdata")
library ('MCMCvis')
library ('coda')
library ('ggplot2')
library('reshape2')
library('bayestestR')
out <- lapply(post, as.mcmc)

# Identify chains with NAs that 
# failed to initialize
NAlist <- c()
for (i in 1:length(out)){
  NAlist[i] <- any (is.na(out[[i]][,1:286]) | out[[i]][,1:286]<0)
}
# Subset chains to those with good initial values
out <- out[!NAlist]
post2 <- post[!NAlist]
outp <- MCMCpstr(out, type="chains")

!NAlist
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE
# default settings for plots 
plt  <- function(object, params,...) {
  MCMCplot(object=out, 
           params=params, 
           guide_axis=TRUE, 
           HPD=TRUE, ci=c(80, 95), horiz=FALSE, 
           #ylim=c(-10,10),
           ...)
  }

Plot model estimates of demographic rates. Life Stages are abbreviated as B = breeder, NB = nonbreeder, FY = first year. First-year abundance accounts for translocated birds.

# Abundance of females at Los Haitises
par(mfrow=c(5,2))
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NFY[",1:13, ", 1]"), 
    main="First-year (FY) minus hacked\n Los Haitises", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,1], 
     ylab="Counts", xlab="Year", type="b")

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NF[",1:13, ", 1]"), 
    main="Adult nonbreeder (NB)\n Los Haitises", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot.new()

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NB[",1:13, ", 1]"),
    main="Adult breeder (B)\n Los Haitises", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot.new()

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NAD[",1:13, ", 1]"), 
    main="Adult Breeders and Nonbreeders\n Los Haitises", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsAdults[,1], 
     ylab="Counts", xlab="Year",  type="b")

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("Ntot[",1:13, ", 1]"), 
    main="All stages\n Los Haitises", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,1]+datl$countsAdults[,1], 
     ylab="Counts", xlab="Year",  type="b")

# Abundance of females at Punta Cana
par(mfrow=c(5,2))
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NFY[",1:13, ", 2]"), 
    main="First-year (FY) plus hacked\n Punta Cana", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,2], 
     ylab="Counts", xlab="Year",  type="b")

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NF[",1:13, ", 2]"), 
    main="Adult nonbreeder (NB)\n Punta Cana", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot.new()

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NB[",1:13, ", 2]"),
    main="Adult breeder (B)\n Punta Cana", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot.new()

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("NAD[",1:13, ", 2]"), 
    main="Adult Breeders and Nonbreeders\n Punta Cana", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsAdults[,2], 
     ylab="Counts", xlab="Year",  type="b")

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("Ntot[",1:13, ", 2]"), 
    main="All stages\n Punta Cana", 
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,2]+datl$countsAdults[,2], 
     ylab="Counts", xlab="Year",  type="b")

Population dynamics are determined by transitions, These transitions between stages are abbreviated as the starting life stage to the final life stage. For example a first-year recruiting to a breeder would be abbreviated as “FY to B”. I’ll list them here for convenience:

“FY to NB” is first-year to nonbreeder.

“NB to NB” is nonbreeder adult to nonbreeder adult.

“B to NB” is a breeding adult to a nonbreeder adult.

“FY to B” is first-year to breeder.

“NB to B” is nonbreeder adult to breeder adult.

“B to B” is breeder adult to breeder adult.

# Finer population segments
par(mfrow=c(4,2))
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 1, ", ", 1:13, ", 1]"), 
    main="Los Haitises\nFirst-years fledged", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 2, ", ", 1:13, ", 1]"), 
    main="Los Haitises\nFY to NB", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 3, ", ", 1:13, ", 1]"), 
    main="Los Haitises\nNB to NB", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 4, ", ", 1:13, ", 1]"), 
    main="Los Haitises\nB to NB", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 5, ", ", 1:13, ", 1]"), 
    main="Los Haitises\nFY to B", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 6, ", ", 1:13, ", 1]"), 
    main="Los Haitises\nNB to B",
    labels = 2011:2023,
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 7, ", ", 1:13, ", 1]"), 
    main="Los Haitises\nB to B", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")


par(mfrow=c(4,2))

plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 1, ", ", 1:13, ", 2]"), 
    main="Punta Cana\nFY fledged",
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 2, ", ", 1:13, ", 2]"), 
    main="Punta Cana\nFY to NB", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 3, ", ", 1:13, ", 2]"), 
    main="Punta Cana\nNB to NB", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 4, ", ", 1:13, ", 2]"), 
    main="Punta Cana\nB to NB", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 5, ", ", 1:13, ", 2]"), 
    main="Punta Cana\nFY to B", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 6, ", ", 1:13, ", 2]"), 
    main="Punta Cana\nNB to B",
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")
plt(object=out, 
    exact=TRUE, ISB=FALSE, 
    params=paste0("N[", 7, ", ", 1:13, ", 2]"), 
    main="Punta Cana\nB to B", 
    labels = 2011:2023, 
    xlab = "Year", ylab= "Abundance")

Other parameter estimates.

# I needed to abbreviate to save plot space
# FY=first-year, NB=Nonbreeder, B=Breeder
par(mfrow=c(1,2))
plt(object=out, 
    params=paste0("mus[",1:8, ", 1]"), 
    exact=TRUE, ISB=FALSE, 
    ylim=c(0,1),
    main="Overall means\n Los Haitises", 
    labels=c("FY survival", "NB survival", "B survival",
             "FY to B", "NB to B", "B to NB",
             "NB detection", "B detection")
    )

plt(object=out, 
    params=paste0("mus[",1:8, ", 2]"), 
    exact=TRUE, ISB=FALSE, 
    ylim=c(0,1),
    main="Overall means\n Punta Cana", 
    labels=c("FY survival", "NB survival", "B survival",
             "FY to B", "NB to B", "B to NB",
             "NB detection", "B detection"))

par(mfrow=c(1,1))
plt(object=out, 
    params="betas", 
    main= "Translocation effects",
    labels=c("FY survival", "NB survival", "B survival",
             "FY to B", "NB to B", "B to NB",
             "NB detection", "B detection"))

# Plot effort effects
plt(object=out, 
    params=paste0("deltas[",1:4,"]"),
    exact=TRUE, ISB=FALSE,
    main= "Effort effects", 
    labels=c("Number of Adults\nEffort", "Number of Adults\nEffort squared",
             "Number of FYs\nEffort", "Number of FYs\nEffort squared"),
    ylim=c(-0.5, 1))

plt(object=out, 
    params=paste0("deltas[",5:8,"]"),
    exact=TRUE, ISB=FALSE,
    main= "Effort effects continued", 
    labels=c("Nonbreeder detection\nEffort", "Nonbreeder detection\nEffort sq",
             "Breeder detection\nEffort", "Breeder detection\nEffort sq"),
    ylim=c(-5, 5))

# Fecundity
par(mfrow=c(1,1))
plt(object=out, 
    params=c("lmu.prod"), 
    labels= c("Productivity\n(log scale)\nLos Haitises",
              "Productivity\n(log scale)\nPunta Cana"))

f <- exp(outp$lmu.prod)

# Is fecundity at LHNP greater than PC
par(mfrow=c(1,1))
fdiff <- f[1,]-f[2,]
hist(fdiff, main="Productivity difference")
abline(v=0, lty=2)

# print probability of direction, similar to frequentist p-value
# so values <=0.025 and >=0.975
# Is the difference in fecundity >0 ?
mean(fdiff>0) 
## [1] 1
# How many times greater is fecundity 
# at treated versus non-treated sites
# median(f.pred[1,]/f[1,])
# median(f.pred[2,]/f[2,])


# gamma = nest treatment effect on fecundity
par(mfrow=c(1,1))
plt(object=out, 
    params=c("gamma"), 
    main="Anti-Parasitic Fly\nTreatment Effects", ylim=c(0,3))

par(mfrow=c(1,1))
sds <- paste0("sds[", 1:9, "]")
plt(object=out, params=sds,
    exact=TRUE, ISB=FALSE,
    main="Temporal SDs (synchrony among sites)",
    labels=c("FY survival", "NB survival", "B survival",
             "FY to B", "NB to B", "B to NB",
             "NB detection", "B detection",
             "Productivity"))

sds2 <- paste0("sds2[", 1:9, "]")
plt(object=out, params=sds2,
    exact=TRUE, ISB=FALSE,
    main="Site-temporal SDs", 
    labels=c("FY survival", "NB survival", "B survival",
             "FY to B", "NB to B", "B to NB",
             "NB detection", "B detection",
             "Productivity"))

# Correlations among vital rates
# Plot is messy with only a few strong correlations
ind <- 1
Rs <- R2s <- c()
for (i in 1:(nrow(outp$R)-1)){
  for (j in (i+1):nrow(outp$R)){
  Rs[ind] <- paste0("R[",i,", ", j, "]")
  R2s[ind] <- paste0("R2[",i,", ", j, "]")
  ind <- ind+1
  }}
par(mfrow=c(2,1))
plt(object=out, params=Rs[1:18], exact=TRUE, ISB=FALSE,
    main="Correlations btw demographic rates\n over time (synchrony)",
    xlab = "Rhos", guide_lines=TRUE)
plt(object=out, params=Rs[19:36], exact=TRUE, ISB=FALSE,
    main="Correlations btw demographic rates\n over time (synchrony), continued...",
    xlab = "Rhos", guide_lines=TRUE)

par(mfrow=c(2,1))
plt(object=out, params=R2s[1:18], exact=TRUE, ISB=FALSE, 
    main="Correlations btw demographic rates\n over time and sites",
    xlab = "Rhos", guide_lines=TRUE)
plt(object=out, params=R2s[19:36], exact=TRUE, ISB=FALSE, 
    main="Correlations btw demographic rates\n over time and sites, continued ...",
    xlab = "Rhos", guide_lines=TRUE)

# Annual averages for integration into the population model
par(mfrow=c(1,1))
labs <- c(paste0("LH ",2011:2023), paste0("PC ",2011:2023))
plt(object=out, params="mn.phiFY", ylim=c(0,1),
    main="First-year survival", labels = labs,
    xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.phiA", ylim=c(0,1),
    main="Adult nonbreeder", labels = labs,
    xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.phiB", ylim=c(0,1),
    main="Breeder", labels = labs,
    xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.psiFYB", ylim=c(0,1),
    main="First-year to breeder", labels = labs,
    xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.psiAB", ylim=c(0,1),
    main="Adult nonbreeder to breeder", labels = labs,
    xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.psiBA", ylim=c(0,1),
    main="Adult breeder to nonbreeder", labels = labs,
    xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.pA", ylim=c(0,1),
    main="Nonbreeder", labels = labs,
    xlab = "Year", ylab= "Detection")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.pB", ylim=c(0,1),
    main="Breeder", labels = labs,
    xlab = "Year", ylab= "Detection")
abline(v=13.5, lwd=2)

plt(object=out, params="mn.prod",
    main="", labels=labs,
    xlab = "Year", ylab= "Productivity")
abline(v=13.5, lwd=2)

3.Model diagnostics

3.1 Check Goodness-of-fit

Goodness-of-fit tests provide evidence that statistical distributions adequately describe the data. Here we test fit for brood size and counts. A Bayesian p-value nearest to 0.5 suggests good fitting statistical distributions, while values near 1 or 0 suggest poor fit.

# Goodness of fit check
fit.check <- function(out, ratio=FALSE, 
                      name.rep="f.dmape.rep", 
                      name.obs="f.dmape.obs",
                      jit=100,
                      ind=1,
                      lab=""){
  par(mfrow=c(1,1))
  # plot mean absolute percentage error
  samps <- MCMCpstr(out, "all", type="chains")
  rep <- samps[name.rep][[1]][ind,]
  obs <- samps[name.obs][[1]][ind,]
  mx <- max(c(rep, obs))
  mn <- min(c(rep, obs))
  plot(jitter(obs, amount=jit), 
       jitter(rep, amount=jit),
       main=paste0("Mean absolute percentage error\n",lab),
       ylab="Discrepancy replicate values",
       xlab="Discrepancy observed values", 
       xlim=c(mn,mx), ylim=c(mn,mx), 
       pch=16, cex=0.5, col="gray10")
  curve(1*x, from=mn, to=mx, add=T, lty=2, lwd=2, col="blue")
  bp1 <- round(mean(rep > obs),2)
  loc <- ifelse(bp1 < 0.5, "topleft", "bottomright")
  legend(loc, legend=bquote(p[B]==.(bp1)), bty="n", cex=2)
  
  if (ratio==TRUE){
    t.rep <- samps["tvm.rep"][[1]][ind,]
    t.obs <- samps["tvm.obs"][[1]][ind,]
    # plot variance/mean ratio
    hist(t.rep, nclass=50,
         xlab="variance/mean ", main=NA, axes=FALSE)
    abline(v=t.obs, col="red")
    axis(1); axis(2)
  }
  return(list('Bayesian p-value'=bp1))
}

# Counts of adults (nonbreeders+breeders)
fit.check(out, ratio=F,
          name.rep="dmape.rep", 
          name.obs="dmape.obs",
          ind=1,
          lab="Adults(Breeder+Nonbreeder)- Poisson", jit=300)

## $`Bayesian p-value`
## [1] 0.67
# Counts of first-year fledglings, ind=2
# Poisson failed fit test bp=0
# So we assigned negative binomial only for Punta Cana.
# Los Haitises excluded from neg bin because no zeroes in counts. 
fit.check(out, ratio=F,
          name.rep="dmape.rep", 
          name.obs="dmape.obs",
          ind=2,
          lab="First-year counts\nNeg binomial-Poisson", jit=300)

## $`Bayesian p-value`
## [1] 0.59
# Productivity
fit.check(out, ratio=F,
          name.rep="f.dmape.rep", 
          name.obs="f.dmape.obs",
          ind=1,
          lab="Productivity-Neg binomial", jit=300)

## $`Bayesian p-value`
## [1] 0.84

3.2 Examine Traceplots, Compare Posterior with Priors

Traceplots provide evidence that models adequately converged.

MCMCtrace(post2, pdf=FALSE, params= "mus", Rhat=TRUE, priors=runif(20000, 0, 1), post_zm=FALSE)

MCMCtrace(post2, pdf=FALSE, params= "betas", Rhat=TRUE, priors=runif(20000, -20, 20), post_zm=FALSE)

MCMCtrace(post2, pdf=FALSE, params= "deltas", Rhat=TRUE, priors=runif(20000, -20, 20), post_zm=FALSE)

MCMCtrace(post2, pdf=FALSE, params= "gamma", Rhat=TRUE, priors=runif(20000, -20, 20), post_zm=FALSE)

MCMCtrace(post2, pdf=FALSE, params= "sds", Rhat=TRUE, priors=rexp(20000, 1), post_zm=FALSE)

MCMCtrace(post2, pdf=FALSE, params= "sds2", Rhat=TRUE, priors=rexp(20000, 1), post_zm=FALSE)

MCMCtrace(post2, pdf=FALSE, params= "R2")